{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 多边形单元上的数值积分 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 问题描述"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "给定一个多边形 $E$ 和 定义在 $E$ 上的函数 $f(\\mathbf x)$，其中 $\\mathbf x = (x, y)\\in E$,  利用数值积分计算:\n",
    "\n",
    "$$\n",
    "\\int_E f(\\mathbf x) \\mathrm d \\mathbf x\n",
    "$$\n",
    "\n",
    "* $f(\\mathbf x)$ 是一个齐次（homogeneous）函数\n",
    "* $E$ 是一个凸或凹的多边形"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 存在方法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Triangulation\n",
    "* Divergence theorem\n",
    "* Moment fitting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 齐次函数\n",
    "\n",
    "$$\n",
    "f(\\lambda\\mathbf x) = \\lambda^qf(\\mathbf x),\\quad \\forall \\lambda > 0\n",
    "$$\n",
    "\n",
    "记 $ \\mathbf y = \\lambda \\mathbf x $, 上式对 $\\lambda$ 求导\n",
    "\n",
    "$$\n",
    "q\\lambda^{q-1}f(\\mathbf x) = \\nabla_{\\mathbf y}f \\cdot\\frac{\\partial\\mathbf y}{\\partial \\lambda} =\\nabla_{\\mathbf y}f \\cdot \\mathbf x\n",
    "$$\n",
    "\n",
    "取 $\\lambda = 1$, 可得：\n",
    "\n",
    "$$\n",
    "qf(\\mathbf x) = \\nabla f(\\mathbf x)\\cdot\\mathbf x\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "给定一个定义在多边形 $E$ 上的齐次函数 $f(\\mathbf x)$, 记 \n",
    "\n",
    "$$\n",
    "\\mathbf F: = \\mathbf x f(\\mathbf x)\n",
    "$$\n",
    "\n",
    "代入散度定理公式\n",
    "\n",
    "$$\n",
    "\\int_E\\nabla\\cdot\\mathbf F \\mathrm d \\mathbf x = \\int_{\\partial E} \\mathbf F\\cdot \\mathbf n\\mathrm ds\n",
    "$$\n",
    "\n",
    "可得\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "& \\int_E\\nabla\\cdot[\\mathbf x f(\\mathbf x)]\\mathrm d \\mathbf x \\\\\n",
    "=& \\int_E (\\nabla\\cdot\\mathbf x)f(\\mathbf x)\\mathrm d \\mathbf x + \n",
    "\\int_E \\mathbf x\\cdot \\nabla f(\\mathbf x)\\mathrm d \\mathbf x\\\\\n",
    "=& 2 \\int_E f(\\mathbf x)\\mathrm d \\mathbf x + \n",
    "\\int_E qf(\\mathbf x)\\mathrm d \\mathbf x\\\\\n",
    "=& (q+2) \\int_E f(\\mathbf x)\\mathrm d \\mathbf x \\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align}\n",
    "& \\int_{\\partial E} (\\mathbf x\\cdot \\mathbf n)  f(\\mathbf x)\\mathrm ds\\\\\n",
    "=& \\sum_{e_i\\in\\partial E}\\int_{e_i} (\\mathbf x\\cdot \\mathbf n_i)  f(\\mathbf x)\\mathrm ds\\\\\n",
    "=& \\sum_{e_i\\in\\partial E}\\int_{e_i} [(\\mathbf x - \\mathbf x_i + \\mathbf x_i)\\cdot \\mathbf n_i] f(\\mathbf x)\\mathrm ds\\\\\n",
    "=& \\sum_{e_i\\in\\partial E}\\int_{e_i} (\\mathbf x_i\\cdot \\mathbf n_i)  f(\\mathbf x)\\mathrm ds\\\\\n",
    "=& \\sum_{e_i\\in\\partial E}(\\mathbf x_i\\cdot \\mathbf n_i)\\int_{e_i} f(\\mathbf x)\\mathrm ds\\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "f(x,y) = x^2 + y^2\n",
    "$$\n",
    "\n",
    "$$\n",
    "f(\\lambda x,\\lambda y) = (\\lambda x)^2 +  (\\lambda y)^2 = \\lambda^2 f(x, y)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gauss 积分"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**代数精度** $\\quad$ 如果求积公式 $\\int_{a}^{b}f(x)\\mathrm dx\\thickapprox \\sum_{k=0}^{n}A_kf(x_k)$, 对所有次数小于等于 $m$ 的多项式是精确的, 但对 $m+1$ 次多项式不精确, 则称其具有 $m$ 次代数精度."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**利用正交多项式构造Gauss求积公式**\n",
    "\n",
    "设 $P_{n}(x),n = 0,1,2\\dots$ 为正交多项式序列, $P_n(x)$ 具有以下性质:\n",
    "* 对每个 $n$, $P_n(x)$ 是 $n$ 次多项式, $n = 0,1,2\\dots$\n",
    "* $\\int_{a}^{b}\\rho(x)p_i(x)p_j(x) \\mathrm dx = 0$\n",
    "* 对任意一个次数小于等于 $n-1$ 的多项式 $P(x)$ ,有 $\\int_{a}^{b}\\rho(x)p(x)P_n(x) \\mathrm dx = 0$, $n\\geq 1$\n",
    "* $P_n(x)$ 在 $(a,b)$ 内有 $n$ 个互异的零点."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**定理** $\\quad$ 设 $x_0,x_1,\\dots,x_n$ 是 $n+1$ 次正交多项式 $P_{n+1}(x)$ 的 $n+1$ 个零点, 则插值型求积公式\n",
    "\n",
    "$$\n",
    "\\int_{a}^{b}\\rho(x)f(x)\\mathrm dx = \\sum_{k=0}^{n}A_kf(x_k)+R_n\n",
    "$$\n",
    "\n",
    "$$\n",
    "A_k = \\int_{a}^{b}\\rho(x)\\Pi_{i=0,i\\neq k}^{n}\\frac{x-x_i}{x_k-x_i}\\mathrm dx\n",
    "$$\n",
    "\n",
    "$$\n",
    "R_n = \\frac{1}{(n+1)!}\\int_{a}^{b}f^{(n+1)}(\\xi(x))\\omega_{n+1}(x)\\rho(x)\\mathrm dx\n",
    "$$\n",
    "\n",
    "是 Gauss 型求积公式.这里 $\\omega_{n+1}(x)$ 是以节点 $x_i(i = 0,1,\\dots,n)$ 为零点的 $n+1$ 次多项式."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**注:**　$n+1$ 次 Legrendre 多项式为\n",
    "\n",
    "$$\n",
    "P_{n+1}(x) = \\frac{1}{(n+1)!2^{n+1}}\\frac{d^{n+1}}{dx^{n+1}}(x^2-1)^{n+1},\\quad x\\in[-1,1],\\quad n = 0,1,2\\dots\n",
    "$$\n",
    "\n",
    "其具有性质:\n",
    "* $n+1$ 次 Legrendre 多项式y与任意不超过 $n$ 次的多项式在区间 $[-1,1]$ 上正交\n",
    "* $n+1$ 次 Legrendre 多项式的 $n+1$ 个零点都在区间 $[-1,1]$ 内"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**利用正交多项式构造 Gauss 求积公式的基本步骤**\n",
    "\n",
    "* 以 $n+1$ 次正交多项式的零点 $x_0,x_1,x_2,\\dots,x_n$ 作为积分点（Gauss点）\n",
    "* 用 Gauss 点对 $f(x)$ 作 Lagrange 插值多项式 $f(x)\\thickapprox \\sum_{k=0}^{n}l_k(x)f(x_k)$，这里 $l_k(x) = \\Pi_{i=0,i\\neq k}^{n}\\frac{x-x_i}{x_k-x_i}$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在区间 $[-1,1]$ 内函数 $f(x)$ 的 Gauss-Legrendre 积分:\n",
    "\n",
    "$$\n",
    "\\int_{-1}^{1}f(x)\\mathrm dx =  \\sum_{i=1}^{n}A_if(x_i)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "其中 $x_i$ 是高斯点, $A_i = \\int_{a}^{b}l_i(x)\\mathrm dx,i = 0,1,2,\\dots,n$ 是权重."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "作变换 $x = \\frac{a+b}{2}+\\frac{b-a}{2}t$ ,即可得到将区间 $[a,b]$ 变换到 $[-1,1]$ 上,所以在任意区间 $[a,b]$ 上函数 $f(x)$ 的 Gauss-Legrendre 积分:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\int_{a}^{b}f(x)\\mathrm dx \n",
    "&= \\frac{b-a}{2} \\int_{-1}^{1}f(\\frac{b-a}{2}x+\\frac{a+b}{2})\\mathrm dx\\\\\n",
    "& = \\frac{b-a}{2}\\sum_{i=1}^{n}A_if(\\frac{b-a}{2}x_i+\\frac{a+b}{2})\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|积分点的个数 $n$ |积分点 $x_i$ | 权重  $A_i$ |\n",
    "|:--------:|:---------:|:--------:|\n",
    "|0  | 0 |2 |\n",
    "|1  | $\\pm 0.5773503$ |1 |\n",
    "|2  | $\\pm 0.7745967$ |0.5555556 |\n",
    "|2  | 0 |0.8888889 |\n",
    "|3  | $\\pm 0.8611363$ |0.3478548 |\n",
    "|3  | $\\pm 0.3399810$ |0.6521452 |\n",
    "|4  | $\\pm 0.9061798$ |0.2369269 |\n",
    "|4  | $\\pm 0.5384693$ |0.4786287 |\n",
    "|4  | 0 |0.5688889 |\n",
    "|5  | $\\pm 0.9324695$ |0.1713245 |\n",
    "|5  | $\\pm 0.6612094$ |0.3607616 |\n",
    "|5  | $\\pm 0.2386192$ |0.4679139 |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gauss–Lobatto rules"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "它与高斯积分相似，具有以下区别:\n",
    "* 积分点包括积分区间的终点\n",
    "* 多项式达到次数是 $2n-3$ 时是精确的, 其中n是积分点的数量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在区间 $[-1,1]$ 内函数 $f(x)$ 的 Lobatto 积分:\n",
    "\n",
    "$$\n",
    "\\int_{-1}^{1}f(x)\\mathrm dx = \\frac{2}{n(n-1)}[f(1)+f(-1)] + \\sum_{i=2}^{n-1}\\omega_if(x_i) + R_n\n",
    "$$\n",
    "\n",
    "其中 $x_i$ 是 $P'_{n-1}(x)$ 的 第 $i-1$ 个零点, $\\omega_i$ 是权重, $R_n$ 是余项."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "权重:\n",
    "\n",
    "$$\n",
    "\\omega_i = \\frac{2}{n(n-1)[P_{n-1}(x_i)]^2},\\qquad x_i \\neq \\pm 1\n",
    "$$\n",
    "\n",
    "这里 $P_{n-1}(x)$ 是 $n-1$ 次 Legrendre 多项式."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "余项:\n",
    "\n",
    "$$\n",
    "R_n = \\frac{-n(n-1)^32^{2n-1}[(n-2)!]^4}{(2n-1)[(2n-2)!]^3}f^{2n-2}(\\xi),\\qquad -1<\\xi<1\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一些权重:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "|积分点的个数 $n$ |$x_i$ | 权重  $\\omega_i$ |\n",
    "|:--------:|:---------:|:--------:|\n",
    "|3  | 0 |$\\frac{4}{3}$ |\n",
    "|3  | $\\pm 1$ |$\\frac{1}{3}$ |\n",
    "|4  | $\\pm\\sqrt{\\frac{1}{5}}$ |$\\frac{5}{6}$ |\n",
    "|4  | $\\pm 1$ |$\\frac{1}{6}$ |\n",
    "|5  | $0$ |$\\frac{32}{45}$ |\n",
    "|5  | $\\pm\\sqrt{\\frac{3}{7}}$ |$\\frac{49}{90}$ |\n",
    "|5  | $\\pm 1$ |$\\frac{1}{10}$ |\n",
    "|6  | $\\pm\\sqrt{\\frac{1}{3}-\\frac{2\\sqrt{7}}{21}}$ |$\\frac{14+\\sqrt{7}}{30}$ |\n",
    "|6  | $\\pm\\sqrt{\\frac{1}{3}+\\frac{2\\sqrt{7}}{21}}$ |$\\frac{14-\\sqrt{7}}{30}$ |\n",
    "|6  | $\\pm 1$ |$\\frac{1}{15}$ |\n",
    "|7  | $0$ |$\\frac{256}{525}$ |\n",
    "|7  | $\\pm\\sqrt{\\frac{5}{11}-\\frac{2}{11}\\sqrt{\\frac{5}{3}}}$ |$\\frac{124+7\\sqrt{15}}{350}$ |\n",
    "|7  | $\\pm\\sqrt{\\frac{5}{11}+\\frac{2}{11}\\sqrt{\\frac{5}{3}}}$ |$\\frac{124-7\\sqrt{15}}{350}$ |\n",
    "|7  | $\\pm 1$ |$\\frac{1}{21}$ |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 数值实验"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "from fealpy.mesh.TriangleMesh import TriangleMesh\n",
    "from fealpy.quadrature.IntervalQuadrature import IntervalQuadrature\n",
    "from fealpy.quadrature.TriangleQuadrature import TriangleQuadrature\n",
    "from meshpy.triangle import MeshInfo, build\n",
    "\n",
    "import matplotlib.pyplot as plt \n",
    "%matplotlib inline\n",
    "from pylab import *\n",
    "\n",
    "def polygon_mesh(points, h):\n",
    "    ## 设置初始网格信息\n",
    "    N = points.shape[0]\n",
    "    facets = np.zeros((N, 2), dtype=np.int)\n",
    "    facets[:, 0] = range(N)\n",
    "    facets[0:-1, 1] = range(1, N)\n",
    "    mesh_info = MeshInfo()\n",
    "    mesh_info.set_points(points)\n",
    "    mesh_info.set_facets(facets.tolist())\n",
    "\n",
    "    ## 生成网格\n",
    "    mesh = build(mesh_info, max_volume=h**2)\n",
    "\n",
    "    ## 转化为 Fealpy 数据结构\n",
    "    point = np.array(mesh.points, dtype=np.float)\n",
    "    cell = np.array(mesh.elements, dtype=np.int)\n",
    "    tmesh = TriangleMesh(point, cell)\n",
    "\n",
    "    return tmesh \n",
    "\n",
    "def exact_quad(f):\n",
    "    tmesh = polygon_mesh(f.points, f.h)\n",
    "    tmesh.uniform_refine(n=4)\n",
    "    \n",
    "    cell = tmesh.ds.cell\n",
    "    point = tmesh.point \n",
    "    \n",
    "    NC = tmesh.number_of_cells()\n",
    "    \n",
    "    qf = TriangleQuadrature(11)\n",
    "    nQuad = qf.get_number_of_quad_points()\n",
    "    b = np.zeros(NC, dtype=np.float)\n",
    "    for i in range(nQuad):\n",
    "        lam_i, w_i = qf.get_gauss_point_and_weight(i)\n",
    "        p = np.einsum('j, ij...->i...', lam_i, point[cell])\n",
    "        b += f(p)*w_i\n",
    "    b *= tmesh.area()\n",
    "    return np.sum(b)\n",
    "    \n",
    "\n",
    "def polygon_quad(f, n):\n",
    "    \n",
    "    qf = IntervalQuadrature(n)\n",
    "    points = f.points\n",
    "    N = points.shape[0]\n",
    "    facets = np.zeros((N, 2), dtype=np.int32)\n",
    "    facets[:, 0] = range(N)\n",
    "    facets[0:-1, 1] = range(1, N)\n",
    "    \n",
    "    norm = points[facets[:, 1]] - points[facets[:, 0]]\n",
    "    W = np.array([[0, -1], [1, 0]])\n",
    "    norm = norm@W\n",
    "    \n",
    "    qf = IntervalQuadrature(n)\n",
    "    nQuad = qf.get_number_of_quad_points()\n",
    "    b = np.zeros(N, dtype=np.float)\n",
    "    for i in range(nQuad):\n",
    "        lam_i, w_i = qf.get_gauss_point_and_weight(i)\n",
    "        p = lam_i[0]*points[facets[:, 0]] + lam_i[1]*points[facets[:, 1]]\n",
    "        b += f(p)*w_i\n",
    "    \n",
    "    b *=np.sum(points*norm, axis=1)\n",
    "    return np.sum(b)/(f.q+2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.2\n",
      "[  1.78267045e-02   7.44949495e-03   8.06818182e-04   1.85528757e-05\n",
      "   1.13797860e-15]\n",
      "[  2.39300847e+00   9.23317684e+00   4.34875000e+01   1.63033608e+10]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEKCAYAAADkYmWmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAACzNJREFUeJzt3XvMZHddx/HPd7t2oaXduqEp27VsLdBV6+0PbdQ/xDZe\ng4BG1EQJpTUk3ogRpWDVQIpiEwUhQU0xpMIftKKpkEZNGgVTuRq1wXAxCC5sKaXttgu9KCDs8Y/z\ne8x0fC5d91mn36evVzLZmXPOnPnNOc++58w5z2ZrmqYA0NOuVQ8AgP87EQdoTMQBGhNxgMZEHKAx\nEQdoTMQfB6rqYFUdr6pt2d9V9XNV9dmqeqCq9lXVg1V14Xas+1Sqqour6vaq+nxV/eKqx3MyquqK\nqvr7VY+D1RPxHaiqDlfV5UuTt+UfBFTV7iSvSfK90zSdPU3T/dM0nTVN0ye3Y/2n2NVJ3jlN095p\nmt6wPLOq3lVVVz2aFVXVDVV17baP8MT4Rx6IOCfsKUn2JPnoqX6hqjptm1d5MMmHt3md2+4UvO+t\nXk8HGrPzdpiqekuSpya5ZZzu+NW1WUmeX1Wfqqp7quqahedUVb28qj5eVfdW1U1Vdc46635Gkn8d\nD49V1d+M6cer6qJxf19V3TJOWXygql619rV/vdM6i0e/4xTBu6vqtVV1NMkrxvSrquojVXVfVf11\nVT11k/f/nKr6UFXdX1XvrKpDY/rfJrksyR+M7fL0LbbjM6vqjqp6SVXdXVV3VtULx7wXJfnpJFeP\ndb1jTN9fVX8+tu8nqurFC+t7QlW9eYzrw1X10qq6Y2H+4aq6uqo+mOShqtpVVS8b++SB8Z5+ZLMx\nL43/bVV1V1Udq6q/q6pvWJh3Q1X9YVX9ZVU9mOR7qur0qvq98fNx15i/Zyx/ztin94x9cEtVnf9o\nx8IpNk2T2w67JTmc5LKFxweTHE9yfZLTk3xzki8kOTTm/1KS9ybZn+SrkvxRkrdusO6DSb6SpBam\nfSXJReP+TUnemvlo/euTHEly29Jzdy08911Jrhr3r0jyX0l+PvMBxp4kz03ysSQXj2nXJHnPBmO7\nOMlDSS5PclqSlyb5tyS7l19rg+cvjuWZYyyvGOv6oSQPJ9k75t+Q5NqF51aSf0zy62P5C5N8PMn3\njfnXjfWfneT8JB9McmRpn/3zmLdnTPuxJOeN+z8+3tt5C9vqtk3eywuTnDH252uT3L4w74Ykx5J8\nx3i8J8nvJ3l7kr1JzkzyjiS/PebvS/KjY7kzk/xpkptX/XPuNvbnqgfgdgp26hyEyxcer8Vz/8K0\nDyT5iXH/I3lk9Pcn+dJibNdZ12KIjye5aET2S0mevjDvVTmxiH9y6fX+KsmVC493jZhesM7YfiPJ\nTQuPK8mnk3z38mttsN2WI/7w0ljvTnLpuL8c8UvXGfvLk7xp3P9E5usIa/N+Zp2IX7HFfr09ybMX\nttWGEV963jljH521MPY/WVrmoSRfu/D4O5P8+wbr+9Yk963659xtvu0Ojyd3L9z/jyRPGvcPJvmL\nqjo+Hlfmo9Dzktx1Aus/N/NR6KcXpt2xwbIbWV7+YJLXV9VrFsY2JTmwzrLnJ/nU2oNpmqZxyuLA\nCY5hzX3TNB1feLy4zZYdTHKgqu5fGOeuJLctjG2r7bI4P1X1giS/nPmoPpmPgp+81aDH6apXJ3ne\nWH4atycneXD59avq3MxH7f9UVWuTd433kKp6YpLXJfmBzB8IleRJVVXTqDqrI+I704n+xTqS+Qj0\nfSf5uvcm+XKSr8l8KiFJLliY//D484zMR37JfKF00fLYjyT5rWmabnwUr/+ZJN+4NO2CLMVxmyyP\n847MR66HNlj+M5m3y9o1hfXO6//POsd5/zdm/ob0vjHt9oywbuGnkjw787exI1W1N/Ppk8XnLo7/\naOYPqEumaVrvQ/tXkjwjybdP03RvVX1L5lM/ax+orJALmzvTZzOf3li02V/+65O8eu2CYVWdW1XP\n2WT5ddc1jlpvTvLKqnpiVX1dkhcszD+a5M7MF1h3jQuaT9vivVyf5Jq1C3NVtbeqnrfBsm9L8qyq\nuqyqdo+Lul9IcrIfTuu5O4/cxv+Q5MFxcfIJVXVaVV1SVd825v9Zkl8bFwkPJPmFLdZ/ZuZTIEfH\ntroy//sDaiNnJfli5ovPZyb5nWwS23E0/cdJXjeOylNVB6rq+xfW959JHqiqfUle+SjHwf8DEd+Z\nrkvym+M3IV4ypi3/JV58/PrMF7JurarPZ77Ieekm699sXS/O/JX7riRvznyR84sL81+U+fe1j2a+\n8Pmezd7INE1vH+/npqr6XJJ/SfKDGyz7sSTPT/KGzN8KnpX5HPKXNxj3Vu9rs/lvSnLJ2MY3jw+w\nH858vvhwknsyh/Hssfy1mT/ADie5NXPUF7fLI157mqaPZv59/Pdn/lC+JMm7txjfmrdk/gZzZ5IP\nZd6fW3lZ5m9P7x/b+dbMF4qT+VTKGZn32XszX6fgMaKc0uJUqqrrMv9GxZWrHstjSVX9bJKfnKbp\nslWPhd4cibOtqupQVX3TuH9p5t/CuHm1o1q9qnpKVX1XzQ5lPs/8uN8unDwXNtluZyW5sar2Zz5v\n/LvTNN2y4jE9Fpye+fz+hUk+l+TGzL+PDyfF6RSAxpxOAWhMxAEaE3GAxkQcoDERB2hMxAEaE3GA\nxkQcoDERB2hMxAEaE3GAxkQcoDERB2hMxAEaE3GAxkQcoDERB2hMxAEaE3GAxnb0f5RclfuTfPWq\nxwHsOMemKftWPYhkh/9HyVWZpim16nEAO8tjqS1OpwA0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMi\nDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4\nQGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIA\njYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0\nJuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCY\niAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMi\nDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4\nQGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIA\njYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0\nJuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCY\niAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMi\nDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4\nQGMiDtCYiAM0JuIAjYk4QGMiDtCYiAM0JuIAjYk4QGMiDtDY7lUP4FSryrTqMQA7zrFVD2BNTZPG\nAXTldApAYyIO0JiIAzQm4gCNiThAYyIO0JiIAzQm4gCNiThAYyIO0JiIAzQm4gCNiThAYyIO0JiI\nAzQm4gCNiThAYyIO0JiIAzQm4gCNiThAYyIO0JiIAzQm4gCNiThAYyIO0JiIAzQm4gCNiThAYyIO\n0JiIAzQm4gCNiThAYyIO0JiIAzQm4gCNiThAYyIO0Nh/A33sMsFpetNUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efbe4d7a128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEPCAYAAABlZDIgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucVXW9//HXm5u3BEUlEwQ11LwTJlLowxEUyTyiloqX\nxGsamqZ20tSO4y+Piae8ZWYeyWNokZdM5Kih4lTeybwmKJ28IWYpIngH/Pz++K6RzWbPzGbYe9be\nM+/n4zGP+e61v2vtz1ow85nv97vW96uIwMzMbFV1yzsAMzPrHJxQzMysIpxQzMysIpxQzMysIpxQ\nzMysIpxQzMysIqqeUCSNkTRb0vOSzijxfi9JUyTNkfSQpIHZ9r6SZkhaJOnyon0OkfSUpCck3SGp\nb7XPw8zMWlfVhCKpG3AFsBewDXCIpM8VVTsGmB8RmwOXAhdl2z8AzgFOLzpm96zebhExBHgaOKlq\nJ2FmZmWpdgtlGDAnIl6KiMXAFGBsUZ2xwHVZ+WZgFEBEvBcRDwIfFtVX9n1tSQJ6A/OqEbyZmZWv\n2gmlP/BKweu52baSdSJiKbCgtS6siFgCTCC1TOYCWwGTKhizmZm1Qy0OyqvVN6UewDeBHSKiPymx\nnNURgZmZWct6VPn4rwIDC14PyLYVmgtsDMzLxkd6R8T8Vo45BIiIeDF7fSOwwmA/gCRPVGZm1g4R\n0eof96VUu4UyExgsaZCkXsA4YGpRnduB8Vn5QGBGieMUntirwNaS1ste7wnMaimAiKj5r3PPPTf3\nGDpDjI7Tcdb6V73E2V5VbaFExFJJJwHTSclrUkTMknQeMDMippHGPyZLmgO8SUo6AEh6AVgb6CVp\nLDA6ImZn+/9J0kfAS8CRLcWwdCl0716lEzQzs09Uu8uLiLgL2LJo27kF5Q+Bg1rYd9MWtl8NXF3O\n519xBZxyStnhmplZO9XioHxF/eAH8Pe/5x1F6xoaGvIOoU31ECM4zkpznJVVL3G2l1alv6zWSYqJ\nE4Pp0+Huu0ErPcRkZtb1SCJqcFA+d6edBgsWwC9+kXckZmadW6dvoUQETz4Je+4JTzwBG22Ud1Rm\nZrXNLZRW7LADnHACTJgAnTh/mpnlqkskFICzz4Y5c+Cmm/KOxMysc+oSXV7NHn4Y9t8fnn4a1l8/\nx8DMzGpYe7u8ulRCATj1VHjjDZg8OaegzMxqnBNKCaUSyrvvwnbbpQce9947p8DMzGqYB+XLtNZa\n8N//nQbpFy7MOxozs86jy7VQmh13HPToAT/7WQcHZWZW49zlVUJrCWXBAth2W7jhBthttw4OzMys\nhrnLayWtsw5ceSUceyy8917e0ZiZ1b8u20JpNm4cDBwIF13UQUGZmdU4d3mVUE5C+ec/011f06bB\nTjt1UGBmZjXMXV7t1K8fXHwxHHMMfPRR3tGYmdWvqicUSWMkzZb0vKQV1n6X1EvSFElzJD0kaWC2\nva+kGZIWSbq8aJ+ekn4u6TlJz0raf1ViPPTQ1O114YWrchQzs66tql1ekroBzwOjgHmkNebHRcTs\ngjrfBLaLiAmSDgb2j4hxktYEhgDbAttGxMkF+zQC3SLiP7LXfSNifonPb7PLq9krr8DQoXDffenu\nLzOzrqpWu7yGAXMi4qWIWAxMAcYW1RkLXJeVbyYlHyLivYh4EPiwxHGPBn7Y/KJUMllZG28M55+f\nur6WLl3Vo5mZdT3VTij9gVcKXs/NtpWsExFLgQWS+rZ0QEl9suL5kh6T9BtJG1Qi2OOOgzXWgMsu\nq8TRzMy6lh55B1BCW82sHsAA4P6IOF3SqcCPgSNKVW5sbPyk3NDQ0Oqazt26wTXXwPDhMHYsfPaz\nKxm5mVkdampqoqmpaZWPU+0xlOFAY0SMyV6fCURETCyoc2dW5xFJ3YHXIqJfwfvjgR2LxlAWRcTa\nWXkAcGdEbFfi88seQyn0ox/BHXfAvfd6HXoz63pqdQxlJjBY0iBJvYBxwNSiOrcD47PygcCMEscp\nPrHbJe2elfcAnq1QvAB8+9vwzjuptWJmZuWp+oONksYAl5GS16SIuFDSecDMiJgmaTVgMvB54E3S\nXWAvZvu+AKwN9AIWAKMjYnZ2a/FkoA/wL+CoiJhb4rPb1UIBeOYZ2H13ePxxGDCgXYcwM6tLflK+\nhFVJKACNjfDYYzB1qru+zKzrqNUur7p21lnw4oswZUrekZiZ1T63UNrw6KOw775pHfoNKnJzsplZ\nbXOXVwmVSCgA3/kOzJsHv/pVBYIyM6txTiglVCqhvPcebL89XHIJ/Nu/VSAwM7Ma5oRSQqUSCqQ5\nvo44It391adP2/XNzOqVE0oJlUwoAMcfn77//OcVO6SZWc1xQimh0gnl7bfTTMS//GV6RsXMrDPy\nbcMdoE8f+NnPvA69mVkpbqG0w2GHwYYbwo9/XPFDm5nlzl1eJVQrofzrX2kd+ttug513rvjhzcxy\n5S6vDrTBBnDppWkxrg9LLf9lZtYFOaG008EHw2abwQ9/2HZdM7OuwF1eq+DVV2HIEJgxI3WBmZl1\nBu7yykH//nDBBXD00bBkSd7RmJnlywllFR17LKy9dhpTMTPrytzlVQH/93/pbq+HH4bBg6v+cWZm\nVVWzXV6SxkiaLel5SWeUeL+XpCmS5kh6KFuNEUl9Jc2QtEjS5S0ce6qkp6p9Dm357Gfh7LNTa+Xj\nj/OOxswsH1VNKJK6AVcAewHbAIdI+lxRtWOA+RGxOXApcFG2/QPgHOD0Fo69P7CwGnG3x8knwwcf\nwH//d96RmJnlo9otlGHAnIh4KSIWA1OAsUV1xgLXZeWbgVEAEfFeRDwIrPCkh6S1gFOB86sV+Mrq\n3h0mTYJzzoFXXsk7GjOzjlfthNIfKPz1OjfbVrJORCwFFkjq28ZxfwD8CHi/QnFWxDbbwLe+Bd/8\nJnTioSkzs5J65B1ACa0OBEnaAfhsRJwmaZO26jc2Nn5SbmhooKGhYZUDbM2ZZ8IXvpBWdzzssKp+\nlJlZRTQ1NdHU1LTKx6nqXV6ShgONETEme30mEBExsaDOnVmdRyR1B16LiH4F748HdoyIk7PXJ5DG\nVj4CegL9gAciYmSJz++Qu7yKzZwJ++yT1qHv16/t+mZmtaRW7/KaCQyWNEhSL2AcMLWozu3A+Kx8\nIDCjxHE+ObGIuCoiBkTEZsAuwHOlkkmedtoJxo9PA/VmZl1FVRNKNiZyEjAd+CswJSJmSTpP0j5Z\ntUnA+pLmAN8GzmzeX9ILwI+B8ZJeLnGHWM067zx47LE0I7GZWVfgBxur6I9/hEMPTevQr7NObmGY\nma0Ur4dSQt4JBWDCBPjoI7jmmlzDMDMrmxNKCbWQUBYuTOvQX3stjBqVayhmZmWp1UH5Lq93b7jq\nKjjuOHj33byjMTOrHrdQOsjXvw7rrw+XXJJ3JGZmrXOXVwm1lFDefDN1ff32t/DFL+YdjZlZy9zl\nVePWWw8uu8zr0JtZ5+WE0oEOPBC23BLOr5kpLc3MKsddXh1s3ry0Dv3dd8MOO+QdjZnZitzlVSc2\n2gguvDB1fXkdejPrTJxQcnDUUbDuunDxxXlHYmZWOe7yyskLL6RJJB98ELbYIu9ozMyWcZdXndl0\nU/j+970OvZl1Hk4oOTrppDSOctVVeUdiZrbq3OWVs1mzYNdd4S9/gYED847GzMxdXnVrq63g1FPh\n+OO9Dr2Z1TcnlBrw3e/Ca6/B9dfnHYmZWftVPaFIGiNptqTnJZ1R4v1ekqZImiPpIUkDs+19Jc2Q\ntEjS5QX115A0TdIsSU9LuqDa51BtPXvCL34B3/kOvP563tGYmbVPVROKpG7AFcBewDbAISWW8T0G\nmB8RmwOXAhdl2z8AzgFOL3Ho/4qIrYDPA7tI2qsa8XekoUPT8ynf+lbekZiZtU+1WyjDgDkR8VJE\nLAamAGOL6owFrsvKNwOjACLivYh4EFhuKsWIeD8i/pCVlwB/AQZU7xQ6zrnnwpNPwq235h2JmdnK\nq3ZC6Q+8UvB6bratZJ2IWAoskNS3nINLWgf4N+DeVQ81f2uskZYKPukkeOutvKMxM1s5PfIOoISy\nblWT1B34FXBpRLzYUr3GxsZPyg0NDTQ0NKxadFW2666w//5w+ulpXMXMrNqamppoampa5eNU9TkU\nScOBxogYk70+E4iImFhQ586sziNZkngtIvoVvD8e2DEiTi469iRgYUSc2srn1/xzKKUsWpQW47rm\nGthzz7yjMbOuplafQ5kJDJY0SFIvYBwwtajO7cD4rHwgMKPEcZY7MUnnA71bSyb1bO214eqr4Rvf\ngHfeyTsaM7PyVP1JeUljgMtIyWtSRFwo6TxgZkRMk7QaMJl0x9abwLjmLixJLwBrA72ABcBoYBFp\nzGUW8BEQwBURsUIHUb22UJodeST06ZNWejQz6yheU76Eek8o8+enrq+bboIRI/KOxsy6ilrt8rJV\n0Lcv/OQnaUbiDz7IOxozs9Y5odS4r34Vtt4afvCDvCMxM2udu7zqwD/+AdtvD7//PXz+83lHY2ad\nnbu8OrENN4SLLoKjj4bFi/OOxsysNCeUOjF+PPTrBz/6Ud6RmJmV5i6vOvLii/CFL8D998PniqfY\nNDOrEHd5dQGbbAKNjV6H3sxqkxNKnZkwIX2/8sp84zAzK+Yurzr03HPpQcc//zm1WszMKsldXl3I\nllum2Yi9Dr2Z1RInlDr1ne/Av/4F113Xdl0zs47gLq869vjjsNde8NRT6VkVM7NK8OSQJXT2hAJw\n9tkwezbcckvekZhZZ+ExlC7q+9+HZ591QjGz/LmF0gk88AAceCA880yaodjMbFXUbAtF0hhJsyU9\nL+mMEu/3kjRF0hxJD0kamG3vK2mGpEWSLi/aZ6ikp7JjXlrtc6h1I0bA174Gp52WdyRm1pVVNaFI\n6gZcAewFbAMcIql40pBjgPkRsTlwKXBRtv0D4Bzg9BKH/hlwTERsAWwhaa9qxF9PLrgA/vAHuOuu\nvCMxs66q2i2UYcCciHgpIhYDU4CxRXXGAs03v94MjAKIiPci4kHgw8LKkjYE1o6ImdmmXwL7VSn+\nuvGpT8HPf56eTVm0KO9ozKwrqnZC6U9a/73Z3GxbyToRsRRYIKm1kYD+2XFaO2aXNHo0jBoF3/te\n3pGYWVdUi3d5rfRAkC3z4x/DrbfCn/6UdyRm1tX0qPLxXwUGFrwekG0rNBfYGJgnqTvQOyLmt3HM\njds45icaGxs/KTc0NNDQ0FBO3HVr3XXhiivSjMRPPAFrrJF3RGZW65qammhqalrl41T1tuEsQTxH\nGhd5DXgUOCQiZhXUmQBsGxETJI0D9ouIcQXvjwe+EBHfKtj2MHAyMBP4X+DyiFhhOLqr3DZcyoEH\nwuDB8MMf5h2JmdWbqj0pnyWFkyPiknYGNga4jNS9NikiLpR0HjAzIqZJWg2YDHweeBMYFxEvZvu+\nAKwN9AIWAKMjYrakHYH/AVYH7oiIU1r47C6bUF5/Pa1Df+edMHRo3tGYWT2p6tQrkh6NiGHtiixH\nXTmhAEyenMZUZs6Enj3zjsbM6kW1E8olQE/gN8C7zdsj4i8r+4EdqasnlAjYe2/YZZc055eZWTmq\nnVDuK7E5ImLkyn5gR+rqCQXg5Zdhxx3hj3+ErbbKOxozqweebbgEJ5Tkyivh+uvTrcTdu+cdjZnV\nuqrO5SWpj6SLJf05+/qxpD4rH6bl4YQToEcP+OlP847EzDqzcru8bgGeYdkUKV8HdoiIA6oY2ypz\nC2WZ55+HL30pDdBvumne0ZhZLav2GMoTETGkrW21xglleRddBHffDdOngzwfgZm1oNrT178vaZeC\nDxsBvL+yH2b5Ou00eOstuPbavCMxs86o3BbKDqRZfZvHTd4CxkfEU1WMbZW5hbKiJ5+EPfdM07Js\ntFHe0ZhZLarmk/LdgK9FxI2SegNExML2hdmxnFBK+/730+qOv/2tu77MbEVV6/KKiI+B72blhfWS\nTKxl55wDzz0HN92UdyRm1pmU2+V1IfAGKz4p39qswLlzC6VlDz0EBxyQWirrrZd3NGZWS6p9l9cL\nJTZHRGy2sh/YkZxQWnfqqfDGG2nOLzOzZtUeQ/liRDzQ3uDy4oTSunffhe22S+un7L133tGYWa2o\ndgvl8Yj4fLsiy5ETStvuvReOOip1ffXunXc0ZlYLqv0cyr2Svir5nqDOZtSotBb9mWfmHYmZ1bty\nWyiLgDWBpcAHpHXfIyJq+m9at1DKs2ABbLst3HAD7LZb3tGYWd6q3ULpAxwJnJ8lkW2APcsMbIyk\n2ZKel3RGifd7SZoiaY6khyQNLHjve9n2WZJGF2w/VdIzkp6SdIOkXmWeh5Wwzjpp4shjj4X3Pf+B\nmbVTuQnlp8Bw4JDs9SLgirZ2ygb0rwD2IiWhQyR9rqjaMcD8iNgcuBS4KNt3a+AgYCvgy8CVSjYC\nvgUMjYjtgR7AOGyVjB2b1k0599y8IzGzelVuQtk5Ik4kdXcREW+R1nlvyzBgTkS8FBGLgSnA2KI6\nY1k2i/HNQPOiXfsCUyJiSbbG/JzseADdgbUk9SB1xc0r8zysFZdfDtddB3/+c96RmFk9KjehLJbU\nHQgASRsAH5exX3/glYLXc7NtJetExFLgbUl9S+z7KtA/IuYBPwZezrYtiIh7yjwPa0W/fnDxxXD0\n0fDRR3lHY2b1ptyEcjlwK9BP0n8C9wMXVCmmVgeCJK1DatUMAjYCPiXp0CrF0uUceihsvDFMnJh3\nJGZWb3qUUykibpD0GDCK9At/v4iYVcaurwIDC14PyLYVmgtsDMzLWkG9I2K+pFez7cX77gH8vXna\nF0m/Bb4E/KpUAI2NjZ+UGxoaaGhoKCPsrkuCq66CoUPT1CzbbJN3RGZWbU1NTTQ1Na3ycaq6pnyW\nIJ4jJaLXgEeBQwqTkaQJwLYRMUHSOFKyGpcNyt8A7Ezq/rob2BzYCZiUff8QuBaYGRErLHDr24bb\n76qr4H/+Bx54wOvQm3U11b5tuF2yMZGTgOnAX0mD7LMknSdpn6zaJGB9SXOAbwNnZvs+C9wIPAvc\nAUyI5FHS4P3jwJOkFtPV1TyPrugb34DVV08D9WZm5ahqCyVvbqGsmr/9DYYPh0cfhc1qehpQM6uk\nmmyhWH0bPDhNyXLcceC8bGZtcUKxVn3727BoEUyalHckZlbr3OVlbXr6aRg5Mi0ZvOuueUdjZtXm\nLi+rmu22gylT4KtfhQrcWWhmnZQTipVl1Ci48UY46KC0hoqZWTEnFCtbQwPccgsccghMn553NGZW\na5xQbKXsuivceiscfjjceWfe0ZhZLXFCsZU2YgRMnQpHHgnTpuUdjZnVCicUa5fhw1MyOeYYuO22\nvKMxs1pQ1uSQZqXstBPccQd85SuwdGmaTNLMui4nFFslO+4Id90FY8akpHLggXlHZGZ5cUKxVTZk\nSLrra6+9YMmSdBeYmXU9TihWEdtvD/fcA6NHp5bK4YfnHZGZdTQnFKuYbbZJSWXPPVNL5cgj847I\nzDqSE4pV1FZbpSfp99gjJZVjj807IjPrKE4oVnFbbgn33ZcmlFy6FI4/Pu+IzKwjVP05FEljJM2W\n9LykM0q830vSFElzJD0kaWDBe9/Lts+SNLpgex9JN2Xb/ypp52qfh62cwYPTRJI//CH8dIXFmc2s\nM6pqC0VSN+AK0pry84CZkm6LiNkF1Y4B5kfE5pIOBi4CmteUPwjYChgA3CNp82w++suAOyLiQEk9\ngDWreR7WPpttlpLK7run7q9TTsk7IjOrpmq3UIYBcyLipYhYDEwBxhbVGQtcl5VvBkZm5X1Ja9Av\niYgXgTnAMEm9gV0j4lqA7P2FVT4Pa6dNNklJ5fLL4eKL847GzKqp2gmlP/BKweu52baSdSJiKfC2\npL4l9n0127Yp8IakayX9RdLVktao1gnYqhs0CP7wB7jqKrjooryjMbNqqcVB+bZWCesBDAVOjIg/\nS7oUOBM4t1TlxsbGT8oNDQ00NDRUJkpbKQMGpJbKyJGweDGcfXbeEZlZs6amJpoqsHpeVZcAljQc\naIyIMdnrM4GIiIkFde7M6jwiqTvwWkT0K64r6S5S0ngReCgiNsu27wKcERH/VuLzvQRwjXnttbRY\n18EHw7kl/wQws7zV6hLAM4HBkgZJ6gWMA6YW1bkdGJ+VDwRmZOWppMH5XpI2BQYDj0bE68ArkrbI\n6o0Cnq3mSVjlfOYz6Zbim26C//gPcL436zyq2uUVEUslnQRMJyWvSRExS9J5wMyImAZMAiZLmgO8\nSUo6RMSzkm4kJYvFwISC5sbJwA2SegJ/B46q5nlYZX360ympjBqV7v76z/8ErfTfQmZWa6ra5ZU3\nd3nVtjfeSNO07LknTJzopGJWK2q1y8usReuvn6ZpufdeOP10d3+Z1TsnFMtV375pQsn7708PPjqp\nmNUvJxTL3brrwt13w8yZcOKJ8PHHeUdkZu3hhGI1oU8f+P3v4amn4IQTnFTM6pETitWM3r3TcsLP\nPZemvV+6NO+IzGxlOKFYTfnUp+COO+DFF+Hoo51UzOqJE4rVnLXWgmnTYN48OOKI9KyKmdU+JxSr\nSWuuCVOnpmdVDj/cScWsHjihWM1aYw247TZYuBAOOSRNKmlmtcsJxWra6qvDrbfCBx+kCSU/+ijv\niMysJU4oVvNWWw1uuSU99Pi1r8GHH+YdkZmV4oRidaFXL7jxxvT9gANSi8XMaosTitWNnj3h17+G\ntdeG/faD99/POyIzK+SEYnWlZ0+4/npYbz3Yd1947728IzKzZk4oVnd69IBf/jIt1rXPPvDuu3lH\nZGbghGJ1qnt3uPZa2GQT2HtveOedvCMys6onFEljJM2W9LykM0q830vSFElzJD0kaWDBe9/Lts+S\nNLpov26S/iKpeElh6yK6d4drroEttoAxY2DRorwjMuvaqppQJHUDrgD2ArYBDpH0uaJqxwDzI2Jz\n4FLgomzfrYGDgK2ALwNXSsut6XcKXku+y+vWDX7+c9h+exg9Gt5+O++IzLquardQhgFzIuKliFgM\nTAHGFtUZC1yXlW8GRmblfYEpEbEkIl4E5mTHQ9IAYG/gmuqGb/WgWzf46U9hp53ScsILFuQdkVnX\nVO2E0h94peD13GxbyToRsRR4W1LfEvu+WrDvJcC/A17fz4C0Hv1ll8GIEbDHHjB/ft4RmXU9PfIO\noAS1+qb0FeCfEfGEpIa26jc2Nn5SbmhooKGhYdUjtJokwcUXw3e/C6NGpaWF11sv76jMal9TUxNN\nTU2rfBxFFRfxljQcaIyIMdnrM4GIiIkFde7M6jwiqTvwWkT0K64r6S7gXFIX2eHAEmANYG3gtxFx\nRInPj2qen9WmCDjrrLSuyj33wAYb5B2RWX2RRES0+sd6KdXu8poJDJY0SFIvYBxQfFfW7cD4rHwg\nMCMrTwXGZXeBbQoMBh6NiLMiYmBEbJYdb0apZGJdlwQXXJAefNx9d3j99bwjMusaqtrlFRFLJZ0E\nTCclr0kRMUvSecDMiJgGTAImS5oDvElKEkTEs5JuJN3JtRiY4OaGlUuCH/wgPQS5++4wYwZsuGHe\nUZl1blXt8sqbu7wM4Pzz03QtM2bARhvlHY1Z7Wtvl1ctDsqbVdQ556Q5wHbbDe67DwYMyDsis87J\nCcW6hDPOSN1fzUll4MC29zGzleOEYl3G6aenpNLQkLq/Ntkk74jMOhcnFOtSTjll+aSy2WZ5R2TW\neTihWJdz4olpYsndd4d774XBg/OOyKxzcEKxLumEE5bdUnzvvWnGYjNbNU4o1mUde2xqqYwcCXff\nDVttlXdEZvXNCcW6tKOOSi2VUaNSUtlmm7wjMqtfTijW5X3966mlssce8Pvfp7VVzGzlOaGYAYce\nmloqe+0Fd94JQ4bkHZFZ/XFCMcscdFBqqYwZk2YqHjo074jM6osTilmBr341tVS+/GWYNi2tAmlm\n5XFCMSsydmxqqeyzD9x2GwwfnndEZvWh2uuhmNWlffaBa69Na6o8+GDe0ZjVBycUsxbsvTdMngz7\n7Qd/+lPe0ZjVPicUs1bstRf86ldwwAFQgSW3zTq1qicUSWMkzZb0vKQzSrzfS9IUSXMkPSRpYMF7\n38u2z5I0Ots2QNIMSX+V9LSkk6t9Dta17bEH3HRTugvs3nvzjsasdlV1xUZJ3YDngVHAPNIa8+Mi\nYnZBnW8C20XEBEkHA/tHxDhJWwM3ADsBA4B7gM2BTwMbRsQTkj4FPAaMLTxmwbG9YqNVzJ/+lO4C\nu/56GD0672jMqqe9KzZWu4UyDJgTES9FxGJgCjC2qM5Y4LqsfDMwMivvC0yJiCUR8SIwBxgWEf+I\niCcAIuIdYBbQv7qnYQa77gq33gqHH54efjSz5VU7ofQHXil4PZcVf/l/UicilgJvS+pbYt9Xi/eV\ntAkwBHikkkGbtWTECJg6FY48Mj2nYmbL1OJzKGU1s7LurpuBU7KWSkmNjY2flBsaGmhoaFjF8Kyr\nGz48JZN99oGrr07PrZjVs6amJpoqcNdJtcdQhgONETEme30mEBExsaDOnVmdRyR1B16LiH7FdSXd\nBZyb1esBTAPujIjLWvl8j6FY1Tz2GHzlK/DTn6axFbPOolbHUGYCgyUNktQLGAdMLapzOzA+Kx8I\nzMjKU4Fx2V1gmwKDgUez934BPNtaMjGrth13hLvuSitA3nhj3tGY5a+qXV4RsVTSScB0UvKaFBGz\nJJ0HzIyIacAkYLKkOcCbpKRDRDwr6UbgWWAxMCEiQtII4DDgaUmPAwGcFRF3VfNczEoZMgSmT0/P\nqyxZkmYtNuuqqtrllTd3eVlHeeaZdCvxxIlpfRWzetbeLq9aHJQ3qzvbbpseetxzT1i6NN0FZtbV\nOKGYVchWW6Wkssceqfvr2GPzjsisYzmhmFXQllvCfffByJGppXL88XlHZNZxnFDMKmzw4DSR5MiR\nqaVy4ol5R2TWMZxQzKpgs81SUtl995RUTjkl74jMqs8JxaxKNtlk+ZbK6afnHZFZdTmhmFXRoEHL\nJ5UzVljAwazzcEIxq7KNN05JZdSolFTOPjvviMyqwwnFrAP075/u/mpOKueem3dEZpXnJ+XNOtDr\nr6ekMmjtz5NNAAAMOklEQVRQem5l4MDlv9ZbD7TSzyebVVZ7n5R3QjHrYG+9lSaVfOUVePnl5b8+\n+GDFJLPxxsuXV1897zOwzs4JpQQnFKs377xTOtE0f82dC+uss2LSKUw8/fpBt2rPI26dmhNKCU4o\n1tl8/HHqNnv55ZYTz8KFMGBA60lnrbXyPhOrZU4oJTihWFf0/vvLkk1LSWfNNUsnnOak85nPQPfu\neZ+J5cUJpQQnFLMVRcAbbyyfYIoTz5tvpqTSWtLp0yfvM7FqqdmEImkMcCnLFtiaWPR+L+CXwI7A\nG8DBEfFy9t73gKOBJaS146eXc8yCYzuhmLXDhx/Cq6+2nHReeim1YFpKOAMHwkYbQc+eeZ+JtUdN\nJhRJ3YDngVHAPNKSwOMiYnZBnW8C20XEBEkHA/tHxDhJWwM3ADsBA4B7gM0BtXXMgmPXRUJpamqi\noaEh7zBaVQ8xguOstJbijIAFC1bsSitMOv/4R7pBoLWks+66lblNut6vZ62p1QW2hgFzIuIlAElT\ngLFA4S//sUDzY143Az/JyvsCUyJiCfBitkTwMFJCaeuYdaUe/pPVQ4zgOCutpTillAzWXRd22KH0\nvkuWwLx5yyecZ59d/pbpjz5qPeEMGACrrdb+OGtNvcTZXtVOKP2BVwpezyUlhZJ1sjXo35bUN9v+\nUEG9V7NtKuOYZpazHj2WJYaWLFy44vjNPfcsK8+bB337Lv8sTvHXBht03DlZ62px6hU/J2zWRfTu\nDdtsk75KWbp02W3SzV9//3uaG625i+2dd1Ldyy9PLafmLrTmcvHrlsodUe/11+F3v6vd+ArL7VHt\nMZThQGNEjMlenwlE4SC6pDuzOo9I6g68FhH9iutKuovUNaa2jllw7NofQDEzq0G1OIYyExgsaRDw\nGjAOOKSozu3AeOAR4EBgRrZ9KnCDpEtIXV2DgUdJd3a1dUygfRfEzMzap6oJJRsTOQmYzrJbfGdJ\nOg+YGRHTgEnA5GzQ/U1SgiAinpV0I/AssBiYkN2yVfKY1TwPMzNrW6d+sNHMzDpO3U8hJ2mSpNcl\nPdVKncslzZH0hKQhHRlfQQytxilpN0kLJP0l+zonhxgHSJoh6a+SnpZ0cgv1cr2e5cRZI9dzNUmP\nSHo8i3OFVVAk9ZI0JbueD0lq5Z6oXOMcL+mfBdfz6I6OM4ujW/b5U0u8l/u1LIiltThr5Vq+KOnJ\n7N/90RbqrNzPekTU9RewCzAEeKqF978M/G9W3hl4uEbj3A2YmvO13BAYkpU/BTwHfK7WrmeZceZ+\nPbM41sy+dwceBoYVvf9N4MqsfDDp2atajHM8cHkNXM9TgetL/dvWyrUsI85auZZ/B9Zt5f2V/lmv\n+xZKRNwPvNVKlbGkqV2IiEeAPpI+3RGxFSojTsj5lumI+EdEPJGV3wFmkW6IKJT79SwzTqiBW9Aj\n4r2suBppzLK4j3kscF1Wvpk0A0SHKyNOyPl6ShoA7A1c00KVmriWZcQJNfB/kxRDazlgpX/W6z6h\nlKH44crmByRr0fCs+fm/2dQzuZG0CalF9UjRWzV1PVuJE2rgemZdH48D/wDujoiZRVWWe7AXWJA9\n2NuhyogT4ICs6+PG7JdmR7sE+HdKJzuokWtJ23FC/tcSUny/lzRT0nEl3l/pn/WukFDqxWPAoIj4\nPHAF8Lu8ApH0KdJfeKdkLYCa1EacNXE9I+LjLIYBwM5lJLZc/nItI86pwCYRMYQ0r951xceoJklf\nAV7PWqaivOvU4deyzDhzvZYFRkTEF0itqRMl7bKqB+wKCeVVYOOC1wOybTUlIt5p7naIiDuBnjn9\npdqD9Et6ckTcVqJKTVzPtuKsletZEM9C4D5gTNFbc8mup9KDvb0jYn4Hh/eJluKMiLciYnH28hrS\n7OAdaQSwr6S/A78Gdpf0y6I6tXAt24yzBq5lcxyvZd//BdzKilNYrfTPemdJKK39xTIVOAI+eXJ/\nQUS83lGBFWkxzsK+SUnDSLd05/GL5RfAsxFxWQvv18r1bDXOWriektaX1CcrrwHsyYqTmDY/2AvL\nP9jbYcqJU9KGBS/Hkp4P6zARcVZEDIyIzUjPqs2IiCOKquV+LcuJM+9rmcWwZtbCR9JawGjgmaJq\nK/2zXotzea0USb8CGoD1JL1Mmp6lF2k6lqsj4g5Je0v6G/AucFQtxgl8TWkq/8XA+6S7VDo6xhHA\nYcDTWX96AGcBg6ih61lOnNTA9QQ+A1yntIxDN+A32fVr88HeGozzZEn7kq7nfODIHOJcQQ1ey5Jq\n8Fp+GrhVaXqqHsANETFd0vGsws+6H2w0M7OK6CxdXmZmljMnFDMzqwgnFDMzqwgnFDMzqwgnFDMz\nqwgnFDMzqwgnFKs4SX2yZ0CaX+8m6fY8Yyqlo+LKHhx8WNJj2TM0he9dLelzbew/tq06lSLphUrM\nKCDpPEkj26izm6QvrupnWe1wQrFqWBeYULStVh94andc2YOA5diDtGzBjhHxwHIfHvGNiCh+er7Y\nfsA27YmxWDYlSWsq8u8UEedGRFtPqjcAX6rE51ltcEKxavghsFm2eNDEbNvakm6SNEvS5OaKkoZK\naspmPL2z1PTYkq6VdJmkByT9TdIB2fblWhiSfiKpeaqIFyRdkM02/Kikz0u6S2mxoG8UHL6PpGmS\nZku6suBYe0p6UNKfJf1G0poFx71Q0p+BrxXFOUjSvUqLFt2ttBDYDsBEYGx2PVYr2uc+SUOz8iJJ\n5yvNQvugpA2yv+D3BS7K9t9U0mbZtZop6Q+Stsj230xpYaknJf1A0qKC6/RHSbcBf8223Zrt/7Sk\nYwtDKvUPmsV2saRnsnNbL9s+JPvMJyTdomVTuFxb8O/0gqTGrIX2pKQtJA0CTgC+nZ3XCElfy+J5\nXFJTqTisxlV60RZ/+Ys0BcpTBa93I60F8xnSL6wHSX+Z9gAeANbL6h0ETCpxvGtJ04EAbAXMKTju\n1IJ6PwGOyMovAN/IyhcDTwBrAusD/yjY/70sXgHTgQOA9YA/AGtk9b4LnFNw3O+0cN5TgcOz8lHA\nrVm5xQWVSBMxDs3KHwN7Z+WJwFkF539AwT73AJ/NysOAe7Py7cBBWfl4YGHBeS4CBhYcY53s++rA\n02QLLWXn17dEnB8D47Ly95vPB3gS2CUrnwdcXBxzdswJWfmbwNVZ+VzgtILPeAr4TFbunff/Y3+t\n/Ffdz+VldePRyGY3lfQEsAnwNrAtcLek5sV+5rWw/+8AImKWpH5lfmZz6+VpYK1Isw+/J+kDSb0L\n4nopi+vXpJU1PwS2Bh7I4upJSoLNftPC530R2D8rTyYlhZXxYUTckZUfI3WVLUdpIr8vATdlsZHF\n1/z5Y7Pyr4D/Ktj10Yh4ueD1tyXtl5UHAJsDJZeBzSwFbszK1wO3ZNewT6TF4yBNw35jqZ1Js9k2\nn9f+LdS5nzSn2I3Ab1uJxWqUE4p1lA8LyktJ//cEPBMRI0rv0uL+zb9Il7B8t+3qLezzcdH+H7Ps\n/37xmEFkx58eEYe1EMu7LWxf1fGHxQXl5mtUrBvwVkQMbePzi7uuPolZ0m7ASGDniPhQ0n2seO3a\n0vxZ5a450nz9WzovImKCpJ2AfYDHJA2NiLZWObUa4jEUq4ZFwNpl1HsO2EBpamwk9VB5Kys2/xJ7\nCdhaUk9J61D+kq+FvwR3zsY+upFmJL6ftKb6CEmfzeJaU9LmZRz3QeCQrHw48Kcy4ykVV6FFQG+A\niFgEvCDpk/EbSdtnxYdZNq7T2ky7fUhJ6UOlu8eGlxFb94JjHwbcH2ntlPladufa10ldheX65Lwg\njQFFxMyIOBf4J8uvxWF1wAnFKi7SuiMPSHpKywbll6uS1VtM+iU1MesGe5zUbVOyfon955K6WJ4B\npgB/aWWflo73KGlFx78C/xcRt0bEG6QpxX8t6UlSotiyjOOeDByVncthwCmt1C0VS0vHngL8ezao\nvWl27GOygfBnSIP2AKcCp2Wf/1lSl2Ipd5EWHPsrcAHwUBkxvAsMk/Q06e6s/5dtHw/8KPvMHQq2\nl3NetwP7Nw/KA/+V/Z95CnggIp5qYT+rUZ6+3qyTkLRGRLyflQ8mDaK3NF6xssdeFBHltDqtC/MY\nilnnsaOkK0hdZ28BR1fw2P7L09rkFoqZmVWEx1DMzKwinFDMzKwinFDMzKwinFDMzKwinFDMzKwi\nnFDMzKwi/j/mLGWOXKt0hwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efbe4d92f28>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "class Example1():\n",
    "    def __init__(self, q=2, h=0.05):\n",
    "        self.h = h\n",
    "        self.q = q\n",
    "        self.points = np.array([\n",
    "                [0.0, 0.0],\n",
    "                [1.0, 0.0],\n",
    "                [1.0, 1.0],\n",
    "                [0.0, 1.0],\n",
    "            ])\n",
    "    def __call__(self, p):\n",
    "        x = p[:, 0]\n",
    "        y = p[:, 1]\n",
    "        q = self.q\n",
    "        return x**q + y**q\n",
    "\n",
    "f = Example1(q=9)\n",
    "fe = exact_quad(f)\n",
    "print(fe)\n",
    "maxit = 5\n",
    "error = np.zeros(maxit, dtype=np.float)\n",
    "fn = np.zeros(maxit, dtype=np.float)\n",
    "for n in range(1, maxit+1):\n",
    "    fn = polygon_quad(f, n)\n",
    "    error[n-1] = np.abs(fn-fe)\n",
    "print(error)\n",
    "print(error[:-1]/error[1:])\n",
    "\n",
    "\n",
    "g = plt.figure()\n",
    "axes = g.gca()\n",
    "axes.set_xlim([-0.1,1.1])  \n",
    "axes.set_ylim([-0.1,1.1])\n",
    "axes.set_axis_off()\n",
    "coord = [[0,0], [1,0], [1,1], [0,1]]\n",
    "coord.append(coord[0])\n",
    "xs, ys = zip(*coord)\n",
    "plt.plot(xs, ys, linewidth=1.0)\n",
    "plt.title('the figure of Integral area')\n",
    "plt.show()\n",
    "\n",
    "x = np.arange(1,maxit+1)\n",
    "y = error\n",
    "plt.plot(x,y)\n",
    "xlabel('the number of integral points')\n",
    "ylabel('error')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.40370198438\n",
      "[  2.61172256e-01   2.82146087e-02   1.95577539e-03   1.18470381e-04\n",
      "   6.65766744e-06   3.54879904e-07   1.81670048e-08   9.00306052e-10\n",
      "   4.34527969e-11   2.07389661e-12]\n",
      "[  9.25663222  14.42630318  16.50855992  17.79457781  18.76033938\n",
      "  19.53431004  20.17869896  20.71917381  20.95224839]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEKCAYAAADkYmWmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAD2hJREFUeJzt3XuwdWVBx/Hv7wVBQcBhJAXCV0l5RbL6o5jUzMDKHPMW\npDNl3honTRknyxtlOmLqTHmbKMfUSGoAxQAH04bJS6gojcXYeCnTuIvcRLmkKL5Pf6x1dLE997P3\nXutZ6/uZOcPee+2z1rPXft/ffs6zfuclpRQkSXXa1fcAJEnbZ4hLUsUMcUmqmCEuSRUzxCWpYoa4\nJFXMEJ+IJLuT7E0yl/c8yQuSfD3JrUkOTXJbkgfOY9+LlOSYJJcl+VaSF/U9np1I8qwkn+h7HOqX\nIT5SSS5PcuLMw3P5pYAk+wJvAn65lHJwKeUbpZSDSilXzGP/C/Yy4KOllENKKafPbkzysSTP3cyO\nkpyR5LVzH+HW+IseE2eIazvuD+wPfGnRB0qyz5x3uRv4wpz3OXcLeN0bHc8sqJRv3AglORN4AHBh\nu9zxRyubgGckuTLJDUlO7XxPkrwiyVeS3JjknCT3WWXfDwH+q717S5J/aR/fm+To9vahSS5slywu\nTXLayo/9qy3rdGe/7RLBJ5O8OclNwKvbx5+b5ItJbk7y4SQPWOf1PynJ55N8I8lHk+xpH/8IcALw\nV+15efAG5/ExSa5O8pIk1ye5Nsmz223PA34beFm7rw+0jx+e5P3t+f1qklM6+7tnkve04/pCkpcm\nubqz/fIkL0vyOeD2JLuSvLx9T25tX9NT1hvzzPjfl+S6JLck+XiSh3W2nZHkr5P8U5LbgF9Ksl+S\nv2j/fFzXbt+/ff592vf0hvY9uDDJEZsdixaolOLXCL+Ay4ETOvd3A3uBdwD7AT8FfAfY025/MXAJ\ncDhwD+DtwFlr7Hs38H0gnce+Dxzd3j4HOItmtn4scBVw8cz37up878eA57a3nwV8D/h9mknG/sCT\ngS8Dx7SPnQp8ao2xHQPcDpwI7AO8FPgfYN/ZY63x/d2xPKYdy6vbfT0euAM4pN1+BvDazvcG+Czw\nx+3zHwh8BfiVdvsb2/0fDBwBfA64auY9+4922/7tYycB92tv/2b72u7XOVcXr/Nang0c0L6fbwYu\n62w7A7gF+Pn2/v7AW4ALgEOAA4EPAH/Wbj8UeGr7vAOB9wLn9f3n3K9iiI/1qw2EEzv3V8Lz8M5j\nlwJPa29/kbuH/uHAd7thu8q+ukG8Fzi6DdnvAg/ubDuNrYX4FTPH+xDwnM79XW2YHrXK2P4EOKdz\nP8A1wC/OHmuN8zYb4nfMjPV64Pj29myIH7/K2F8BvLu9/VWa6wgr2353lRB/1gbv62XAEzvnas0Q\nn/m++7Tv0UGdsf/dzHNuBx7Uuf8I4H/X2N/PADf3/efcr8K+aGqu79z+P+De7e3dwPlJ9rb3QzML\nvR9w3Rb2fxjNLPSazmNXr/Hctcw+fzfwtiRv6oytAEeu8twjgCtX7pRSSrtkceQWx7Di5lLK3s79\n7jmbtRs4Msk3OuPcBVzcGdtG56W7nSTPBP6AZlYPzSz4vhsNul2uej1wcvv80n7dF7ht9vhJDqOZ\ntf97kpWHd7WvgST3At4KPI7mAyHAvZOktKmufhji47XVv1hX0cxAP73D494I3AX8OM1SAsBRne13\ntP89gGbmB82F0q7ZsV8FvK6UcvYmjv814CdnHjuKmXCck9lxXk0zc92zxvO/RnNeVq4prLau/4N9\ntuv+f0PzE9Kn28cuow3WDfwW8ESan8auSnIIzfJJ93u747+J5gPquFLKah/afwg8BPi5UsqNSX6a\nZuln5QNVPfHC5nh9nWZ5o2u9v/zvAF6/csEwyWFJnrTO81fdVztrPQ94TZJ7JXko8MzO9puAa2ku\nsO5qL2j+xAav5R3AqSsX5pIckuTkNZ77PuAJSU5Ism97Ufc7wE4/nFZzPXc/x/8G3NZenLxnkn2S\nHJfkZ9vt5wKvbC8SHgm8cIP9H0izBHJTe66ew49+QK3lIOBOmovPBwJvYJ2wbWfT7wTe2s7KSXJk\nkl/t7O/bwK1JDgVes8lxaMEM8fF6I/Cqtgnxkvax2b/E3ftvo7mQdVGSb9Fc5Dx+nf2vt69TaH7k\nvg54D81Fzjs7259H09e+iebC56fWeyGllAva13NOkm8C/wn82hrP/TLwDOB0mp8KnkCzhnzXGuPe\n6HWtt/3dwHHtOT6v/QD7dZr14suBG2iC8eD2+a+l+QC7HLiIJtS75+Vuxy6lfImmj/8Zmg/l44BP\nbjC+FWfS/ARzLfB5mvdzIy+n+enpM+15vojmQjE0SykH0Lxnl9Bcp9AAxOUsLVqSN9I0Kp7T91iG\nJMnzgaeXUk7oeyyqlzNxzV2SPUke3t4+nqaFcV6/o+pfkvsneWQae2jWmSd/XrQzXtjUIhwEnJ3k\ncJp14z8vpVzY85iGYD+a9f0HAt8Ezqbp40vb5nKKJFXM5RRJqpghLkkVM8QlqWKGuCRVzBCXpIoZ\n4pJUMUNckipmiEtSxQxxSaqYIS5JFTPEJalihrgkVcwQl6SKGeKSVDFDXJIqZohLUsUMcUmqmCEu\nSRXz/7GpH0g4BtjT9zikyl1WCtcs62CGuABIeDLwLuBSwP/xqrQ9DwI+Arx4WQc0xEXC84E/BR5f\nCp/tezxSrRJOAY5Z5jEN8QlLCHAa8HTg0aXw1Z6HJGmLDPGJSrgH8E7gWOCRpXBjz0OStA2G+AQl\nHAScC3wPOLEU7uh5SJK2yYrhxCTcH/g4cBXwVANcqpshPiEJe4BLgAuA3yuFu3oekqQdcjllIhIe\nAZwPnFoKf9v3eCTNhyE+AZ0O+DNL4cN9j0fS/LicMnJtB/ztNB1wA1waGWfiI2UHXJoGQ3yE7IBL\n02GIj4wdcGlaXBMfETvg0vQY4iNhB1yaJpdTRsAOuDRdhnjl7IBL0+ZySsXsgEtyJl4hO+CSVhji\nlbEDLqnLEK+IHXBJs1wTr4QdcEmrMcQrYAdc0lpcThk4O+CS1mOID5gdcEkbcTlloOyAS9oMZ+ID\nYwdc0lYY4gNiB1zSVhniA2EHXNJ2uCY+AHbAJW2XId4zO+CSdsLllB7ZAZe0U4Z4T+yAS5oHl1N6\nYAdc0rw4E18iO+CS5s0QXxI74JIWwRBfAjvgkhbFNfEFswMuaZEM8QWyAy5p0VxOWRA74JKWwRBf\nADvgkpbF5ZQ5swMuaZmcic+JHXBJfTDE58AOuKS+GOI7ZAdcUp9cE98BO+CS+maIb5MdcElD4HLK\nNtgBlzQUhvgW2QGXNCQup2yBHXBJQ+NMfBPsgEsaKkN8A3bAJQ2ZIb6OtgP+fuC72AGXNECuia+h\n7YD/K3AldsAlDZQhvopOB/x87IBLGjCXU2bYAZdUE0O8I+EpNBcx7YBLqoIh3kp4AfAqmg74Z/se\njyRtxuRDvO2Avw54GnbAJVVm0iHedsDfBTwUO+CSKjTZELcDLmkMJlkxtAMuaSwmF+J2wCWNyaSW\nU+yASxqbyYS4HXBJYzSJELcDLmmsRh3idsAljd1oQ9wOuKQpGGWI2wGXNBWjqxjaAZc0JaMKcTvg\nkqZmNMspdsAlTdEoQtwOuKSpqj7E7YBLmrJqQ9wOuCRVGuJ2wCWpUV2I2wGXpB+qqmJoB1yS7q6a\nELcDLkk/qorlFDvgkrS6wYe4HXBJWtugQ9wOuCStb5AhbgdckjZncCFuB1ySNm9QIW4HXJK2ZjAV\nQzvgkrR1gwhxO+CStD29L6fYAZek7es1xO2AS9LO9BbidsAlaeeWHuJ2wCVpfpYa4nbAJWm+lj0T\n/yB2wCVpbpZdMXwE8A8GuCTNx7Jn4o8BPpjwY6Xwl0s+tiSNzlJDvBQuS3gU8M8JRwGvKIW9yxyD\nJI3J0n9jsxSuAB4F/ALw9wn7LXsMkjQWvfzafSncDDwWOAD4UMIhfYxDkmrX27+dUgrfBk4G/hu4\nOOGIvsYiSbXq9R/AKoXvAy8CzgEuSXhYn+ORpNr0/g9glUIB3pBwLfCxhJNL4RN9j0uSajCIf4oW\noBTOBH4H+MeEk/oejyTVoPeZeFcpXJTwOJou+RF2ySVpfYMKcbBLLklbMZjllC675JK0OYMMcbBL\nLkmbMdgQB7vkkrSRQYc42CWXpPUM7sLmauySS9LqBj8T77JLLkl3V8VMvMsuuST9UHUhDnbJJWlF\nVcspXXbJJaniEAe75JJUdYiDXXJJ01Z9iINdcknTVeWFzdXYJZc0RaOYiXfZJZc0JaOZiXfZJZc0\nFaMMcbBLLmkaRrec0mWXXNLYjTrEwS65pHEbfYiDXXJJ4zWJEAe75JLGabQXNldjl1zS2ExmJt5l\nl1zSWExqJt5ll1zSGEw2xMEuuaT6TXI5pcsuuaSaTT7EwS65pHoZ4i275JJqZIh32CWXVJtJX9hc\njV1ySTVxJr4Gu+SSauBMfB12ySUNnSG+AbvkkobM5ZRNsEsuaagM8U2ySy5piAzxLbBLLmloDPEt\nsksuaUi8sLkNdsklDYUz8R2wSy6pb87Ed8guuaQ+GeJzYJdcUl9cTpkTu+SS+mCIz5FdcknLZojP\nmV1ySctkiC+AXXJJy+KFzQWxSy5pGZyJL5hdckmL5Ex8CeySS1oUQ3xJ7JJLWgSXU5bILrmkeTPE\nl8wuuaR5MsR7YJdc0rwY4j2xSy5pHryw2SO75JJ2ypn4ANgll7RdzsQHwi65pO0wxAfELrmkrXI5\nZWDskkvaCkN8gOySS9osQ3yg7JJL2gxDfMDskkvaiBc2B84uuaT1OBOvhF1ySatxJl4Ru+SSZhni\nlbFLLqnL5ZQK2SWXtMIQr5RdcklgiFfNLrkkQ7xyq3TJj+15SJKWyAubIzDTJf94wkml8Mm+xyVp\n8ZyJj0inS35ewm/0PR5Ji+dMfGRW6ZKf3veYJC2OIT5Cq3TJX2mXXBonl1NGqtMlfzRwpl1yaZwM\n8RHrdMkPpOmSH9zzkCTNmSE+cnbJpXEzxCeg0yV/L3bJpVHxwuZE2CWXxsmZ+MTYJZfGxZn4BNkl\nl8bDEJ8ou+TSOLicMmF2yaX6OROfuFK4OeGxwFnADQkvBG7teVhSrR4O3LnMA6aUsszjaaAS9gHu\nAm4BPtXzcKSanVkK5y7rYIa4JFXMNXFJqpghLkkVM8QlqWKGuCRVzBCXpIoZ4pJUMUNckipmiEtS\nxQxxSaqYIS5JFTPEJalihrgkVcwQl6SKGeKSVDFDXJIqZohLUsUMcUmqmCEuSRUzxCWpYoa4JFXM\nEJekihniklQxQ1ySKmaIS1LFDHFJqpghLkkVM8QlqWKGuCRVzBCXpIoZ4pJUMUNckipmiEtSxQxx\nSaqYIS5JFTPEJalihrgkVcwQl6SKGeKSVLH/B9dew3HC8/mPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efc0e644940>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEPCAYAAABRHfM8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH5tJREFUeJzt3XucVWXd9/HPl4MKKoiCmiKk4Qm988xYWk6phT0+YdaT\nWJadD2YHe3U/aq+S4bbDrZ3Luq1X3t49PhZiajd2BA9jTx4AAUWRU4oIgmYCSqKGw+/5Y10jm2GG\nWXNYs/bs/X2/XvOatde+1lq/PTD7N+v67eu6FBGYmZl1ZkDZAZiZWf/ghGFmZrk4YZiZWS5OGGZm\nlosThpmZ5eKEYWZmuRSeMCRNlLRE0jJJF7fz/CclLZS0QNKfJR1W8dylkpZLWizpbUXHamZmHVOR\n4zAkDQCWAacCa4C5wOSIWFLRZreI+Efa/p/ABRFxhqTxwPXACcBo4Dbg4PDAETOzUhR9hzEBWB4R\nKyNiMzANmFTZoDVZJLsBW9L2O4FpEfFKRDwOLE/nMzOzEgwq+Pz7A6sqHq+mnTd9SRcAXwQGA2+t\nOPbeimZPpn1mZlaCqih6R8RPImIccDHw1bLjMTOz7RV9h/EkMKbi8ei0ryM3AFdXHHtAZ8dKck3D\nzKwbIkJdaV/0HcZcYJyksZJ2AiYDMyobSBpX8fBMsiI5qd1kSTtJOhAYB8xp7yIRUXVfU6ZMKT0G\nx+SY6jEux5TvqzsKvcOIiBZJFwIzyZLTNRGxWNJUYG5E/Ba4UNJpwD+B9cD56dhHJE0HHgE2k316\nyncTZmYlKbpLioj4I3Bom31TKra/sINjvwl8s7jozMwsr6ooeteixsbGskPYjmPKxzHlV41xOabi\nFDpwry9Ick+VmVkXSSKqrOhtZmY1wgnDzMxyccIwM7NcnDDMzCwXJwwzM8vFCcPMzHJxwjAzs1yc\nMMzMLBcnDDMzy8UJw8zMcnHCMDOzXJwwzMwsFycMMzPLxQnDzMxyccIwM7NcnDDMzCwXJwwzM8vF\nCcPMzHJxwjAzs1ycMMzMLBcnDDMzy8UJw8zMcqmJhNHSUnYEZma1ryYSxiOPlB2BmVntq4mEMXt2\n2RGYmdU+JwwzM8ul8IQhaaKkJZKWSbq4necvkrRI0gOSZkk6oOK5FknzJS2Q9JuOruGEYWZWPEVE\ncSeXBgDLgFOBNcBcYHJELKlocwowOyJekvQpoDEiJqfnno+IYZ1cI3bdNXjqKdhtt8JeiplZTZFE\nRKgrxxR9hzEBWB4RKyNiMzANmFTZICLuioiX0sP7gP0rns71Yl7/erj//t4I18zMOlJ0wtgfWFXx\neDXbJoS2Pgr8oeLxzpLmSLpH0qSODmpocLeUmVnRBpUdQCtJ5wHHAadU7B4bEWslHQjcIWlhRKxo\ne2xDA0yf3leRmpnVp6ITxpPAmIrHo9O+bUg6DbgUeHPqugIgItam7yskNQPHANsljHvvbWLWLGhq\ngsbGRhobG3vzNZiZ9XvNzc00Nzf36BxFF70HAkvJit5rgTnAuRGxuKLNMcCNwNsj4tGK/XsAmyLi\nn5JGAncDkyoL5qldbNkS7L03LFgAo0cX9nLMzGpG1RW9I6IFuBCYCSwCpkXEYklTJZ2Zml0J7Arc\n2Objs4cD90taANwOfLNtsmgluY5hZla0Qu8w+oKkiAguvxw2boQrryw7IjOz6ld1dxh9yXcYZmbF\nqpk7jA0b4IADYP16GFQ1n/0yM6tOdX2HsccesP/+nrnWzKwoNZMwwN1SZmZFcsIwM7NcnDDMzCyX\nmil6A2zeDCNGwNq1sPvuJQdmZlbF6rroDTB4MBx1lGeuNTMrQk0lDHC3lJlZUZwwzMwsl5pNGP28\nNGNmVnVqLmGMHQstLbB6ddmRmJnVlppLGJ651sysGDWXMMAJw8ysCDWZMCZMcMIwM+ttNTVwr9WG\nDdnKexs2eOZaM7P21P3AvVZ77JFNdf7ww2VHYmZWO2oyYYDrGGZmvc0Jw8zMcnHCMDOzXGqy6A1b\nZ65dswaGDSshMDOzKuaid4XBg+Hoo2Hu3LIjMTOrDTWbMCDrlpozp+wozMxqQ80nDNcxzMx6R10k\njH5epjEzqwo1nTDGjMmSxapVZUdiZtb/1XTC8My1Zma9p6YTBjhhmJn1FicMMzPLpfCEIWmipCWS\nlkm6uJ3nL5K0SNIDkmZJOqDiufPTcUslfbA71z/hBFiwIBvIZ2Zm3VdowpA0ALgKeDtwBHCupMPa\nNJsPHBcRRwM3Ad9Kx44ALgNOABqAKZKGdzWGYcOyZVs9c62ZWc8UfYcxAVgeESsjYjMwDZhU2SAi\n7oqIl9LD+4D90/bbgZkR8VxEbABmAhO7E4S7pczMeq7ohLE/UPmh1tVsTQjt+Sjwhw6OfbKTYzvk\nFfjMzHquatajk3QecBxwSlePbWpqenW7sbGRxsbGbZ5vaIAf/rBn8ZmZ9WfNzc00Nzf36ByFzlYr\n6USgKSImpseXABERV7RpdxrwA+DNEfFs2jcZaIyIT6XHVwN3RsQNbY5td7baSq+8kq3C9+STMLzL\nVRAzs9pTjbPVzgXGSRoraSdgMjCjsoGkY4CrgXe2JovkT8DpkoanAvjpaV+XDRoExxzjmWvNzHqi\n0IQRES3AhWQF60XAtIhYLGmqpDNTsyuBXYEbJS2Q9Jt07HrgcuB+YDYwNRW/u8WFbzOznqnZBZTa\nuvFGuO46mDGj06ZmZjWvGrukqoZnrjUz65m6SRgHHAADBsDKlWVHYmbWP9VNwvDMtWZmPVM3CQOc\nMMzMeqLuEobX+DYz6566+ZQUwPPPw377wfr1MHhwwYGZmVUxf0qqE8OGwWtfCw89VHYkZmb9T10l\nDHAdw8ysu5wwzMwsFycMMzPLpa6K3pDNXDtiBKxalc1ga2ZWj1z0zmHQIDj2WM9ca2bWVXWXMMAr\n8JmZdUddJgzXMczMuq6uE0Y/L9+YmfWpukwYo0dntYzHHy87EjOz/qMuE4ZnrjUz67q6TBjghGFm\n1lVOGGZmlkvdDdxrtXEj7LtvNnPtTjsVEJiZWRXzwL0u2H13OOggWLiw7EjMzPqHuk0Y4G4pM7Ou\n6DRhSBoo6aK+CKavOWGYmeXXacKIiBbg3D6Ipc95yVYzs/xyFb0lfQ8YDNwAvNC6PyLmFxdaPt0t\negO0tGQz165cmX03M6sX3Sl6D8rZ7uj0/d8q9gXw1q5crNoMHLh15tq3va3saMzMqluuhBERbyk6\nkLK01jGcMMzMdizXp6QkDZf0XUn3p6/vSBpedHB9wYVvM7N88n6s9j+BjcB709fzwLV5DpQ0UdIS\nScskXdzO82+SNE/SZklnt3muRdJ8SQsk/SZnrF3imWvNzPLJW/R+ICKO7mxfO8cNAJYBpwJrgLnA\n5IhYUtFmDDAM+BIwIyJurnju+YgY1sk1ul30bjV6NPz5z9lAPjOzelDkSO8XJZ1ccaGTgBdzHDcB\nWB4RKyNiMzANmFTZICKeiIiHyYrobXXpxXSXV+AzM+tc3oTxKeDHkh6X9DhwFfDJHMftD6yqeLw6\n7ctrZ0lzJN0jaVLnzbvHdQwzs851+imp1K10aEQcJWkYQEQ8X3hkmbERsVbSgcAdkhZGxIq2jZqa\nml7dbmxspLGxsUsXaWiASy/tYaRmZlWsubmZ5ubmHp0jbw3j/og4vssnl04EmiJiYnp8CRARcUU7\nba8Fbq2sYeR5vjdqGP/4B+yzj2euNbP6UWQN4zZJX5J0gKQ9W79yHDcXGCdprKSdgMnAjB20fzV4\nSXukY5A0Engj8EjOeLtkt93gda+DBx8s4uxmZrUh70jvc9L3z1TsC2CHnyuKiBZJFwIzyZLTNRGx\nWNJUYG5E/FbS8cAtwB7AmZKaIuJfgMOBn0pqScd+s/LTVb2ttY5xwglFXcHMrH/rtEsq1TDeEBF3\n901IXdMbXVIAP/853HUXXHddLwRlZlblCumSiogtZJ+Kqmn+pJSZ2Y7lrWHcLundkvpkXEQZxo+H\np56CdevKjsTMrDrlTRifBKYDL0t6XtJGSX310do+MXAgHHec18cwM+tI3oQxHPgQ8LU0VccRwOlF\nBVUWd0uZmXUsb8L4MXAiW1fe20gN1jWcMMzMOpY3YTRExGeAlwAiYj1Qc0PcWpds9cy1Zmbby5sw\nNksaSJogUNIoYEthUZVkv/1gyBB47LGyIzEzqz55E8YPyQbX7S3p68BfgG8UFlWJ3C1lZta+vEu0\nXi9pHtm6FgLOiojFhUZWktaE8b73lR2JmVl1yTs1CGlajsKm5qgWDQ1w001lR2FmVn1yzVZbzXpr\napBWrTPXrlsHO+/ca6c1M6sqRc5WWzd22w3GjfPMtWZmbTlhtMOFbzOz7TlhtMMJw8xse04Y7XDC\nMDPbnove7WhpgREjYMUK2GuvXj21mVlVcNG7lwwcCMcf75lrzcwqOWF0wN1SZmbbcsLogBOGmdm2\nXMPowNq1cOSR8Pe/Q+2uM2hm9co1jF70mtfArrvCX/9adiRmZtXBCWMH3C1lZraVE8YOOGGYmW3l\nhLEDThhmZlu56L0DmzbBqFHw7LOwyy6FXMLMrBQueveyoUPhkEPggQfKjsTMrHxOGJ1oaPCIbzMz\ncMLolOsYZmaZwhOGpImSlkhaJunidp5/k6R5kjZLOrvNc+en45ZK+mDRsbZnwgQnDDMzKLjoLWkA\nsAw4FVgDzAUmp/XBW9uMAYYBXwJmRMTNaf8I4H7gWEDAPODYiHiuzTUKK3pDNnPtnnvCo4/CyJGF\nXcbMrE9VY9F7ArA8IlZGxGZgGjCpskFEPBERDwNt3/XfDsyMiOciYgMwE5hYcLzb8cy1ZmaZohPG\n/sCqiser077uHPtkF47tVa5jmJnBoLID6A1NTU2vbjc2NtLY2Nir529ogP/4j149pZlZn2pubqa5\nublH5yi6hnEi0BQRE9PjS4CIiCvaaXstcGtFDWMy0BgRn0qPrwbujIgb2hxXaA0D4KmnYPz4bACf\nZ641s1pQjTWMucA4SWMl7QRMBmbsoH1l8H8CTpc0PBXAT0/7+ty++8Luu8Py5WVc3cysOhSaMCKi\nBbiQrGC9CJgWEYslTZV0JoCk4yWtAt4DXC3poXTseuBysk9KzQampuJ3KVzHMLN657mkcvrOd2DF\nCrjqqsIvZWZWuGrskqoZvsMws3rnO4ycNm3KBu6tW+eZa82s//MdRoGGDoXDDoMFC8qOxMysHE4Y\nXeBuKTOrZ04YXeCEYWb1zAmjC5wwzKyeOWF0waGHZkXvZ54pOxIzs77nhNEFAwbACSf4LsPM6pMT\nRhe5W8rM6pUTRhdNmOC1McysPnngXhc9/XQ2HuPZZ7MuKjOz/sgD9/rAPvvA8OGeudbM6o8TRje4\njmFm9cgJoxucMMysHjlhdIMThpnVIxe9u+HFF2GvvbLC95AhfXppM7Ne4aJ3HxkyBA4/3DPXmll9\nccLoJndLmVm9ccLoJicMM6s3Thjd5IRhZvXGCaObDjkENmyAv/2t7EjMzPqGE0Y3eeZaM6s3Thg9\n4G4pM6snThg94IRhZvXEA/d64Jln4OCDs1X4PHOtmfUnHrjXx0aNgj33hKVLy47EzKx4Thg95G4p\nM6sXThg9NGGCE4aZ1YfCE4akiZKWSFom6eJ2nt9J0jRJyyXdK2lM2j9W0iZJ89PXT4qOtTsaGrxk\nq5nVh0KL3pIGAMuAU4E1wFxgckQsqWjzaeBfIuICSecA74qIyZLGArdGxOs7uUZpRW/wzLVm1j9V\nY9F7ArA8IlZGxGZgGjCpTZtJwC/S9q/JkkurLr2YMgwZAuPHw/z5ZUdiZlasohPG/sCqiser0752\n20REC7BB0p7puddKmifpTkknFxxrt7nwbWb1YFDZAbSj9a5iLTAmItZLOhb4jaTxEfGPEmNrV0MD\n/O53ZUdhZlasohPGk8CYisej075Kq4EDgDWSBgLDImJdeu6fABExX9KjwCHAdp0/TU1Nr243NjbS\n2NjYS+Hn09AAl13Wp5c0M+uS5uZmmpube3SOooveA4GlZHWJtcAc4NyIWFzR5gLgyFT0ngyclYre\nI4F1EbFF0kHAXWTF8Q1trlFq0Rtgy5as8L1kCeyzT6mhmJnlUnVF71STuBCYCSwCpkXEYklTJZ2Z\nml0DjJS0HPgCcEna/2ZgoaT5wHTgk22TRbUYMMDjMcys9nkuqV5y2WXQ0gJf/3rZkZiZda7q7jDq\niT8pZWa1zncYveSZZ2DcOFi/3jPXmln18x1GiUaNgpEjs8K3mVktcsLoRe6WMrNa5oTRi5wwzKyW\nOWH0IicMM6tlLnr3opdeygbwPfMMDB1adjRmZh1z0btku+wCRxwB8+aVHYmZWe9zwuhlHvFtZrXK\nCaOXuY5hZrXKCaOXOWGYWa1ywuhlBx8ML7wAN94IVVKLNzPrFU4YvUyC6dOzSQgbGuC228qOyMys\nd/hjtQXZsiVLHF/9KowZA9/4RpZAzMyqgT9WW0UGDIDJk+GRR7Lv73kPnHUWPPxw2ZGZmXWPE0bB\nBg+Gj38cli+HN78ZTj0VPvABeOyxsiMzM+saJ4w+sssu8MUvZolj3LhsvMYFF8CaNWVHZmaWjxNG\nHxs2DKZMyaZBHzoUjjwSLr4Ynn227MjMzHbMCaMkI0fCt78NCxfChg1w6KHwta/Bxo1lR2Zm1j4n\njJKNHg0//Sncdx8sXpyN4/j+97OJDM3MqokTRpUYNw6uvx5mzoTbb4dDDoFrroFXXik7MjOzjMdh\nVKl77oEvfxnWroXLL88+luu1ws2st3RnHIYTRhWLgFmzssSxZUs2enzixGw0uZlZTzhh1KgIuPlm\n+MpXYNSobNT4ySeXHZWZ9Wce6V2jJHj3u+Ghh+AjH4HzzoN3vAMWLCg7MjOrJ04Y/cigQfChD8HS\npXDGGVnSOOccWLas7MjMrB44YfRDO+8Mn/0s/PWvcNRR8MY3wsc+BqtWlR2ZmdUyJ4x+bNdds4L4\n8uWw995w9NFw0UXwzDNlR2ZmtajwhCFpoqQlkpZJurid53eSNE3Sckn3ShpT8dylaf9iSW8rOtb+\nasSIrBC+aFE2buOww+Cyy+C558qOzMxqSaEJQ9IA4Crg7cARwLmSDmvT7KPAuog4GPg+cGU6djzw\nXuBw4AzgJ1L/+UBpc3Nzn19z333hRz+C+++HJ57IRo1/61vw4ovlxdQZx5RPNcYE1RmXYypO0XcY\nE4DlEbEyIjYD04BJbdpMAn6Rtn8NvDVtvxOYFhGvRMTjwPJ0vn6hzP8gBx4I//Vf0NycTTkybhxc\nfTXcfHMza9dmc1e9/HJ1LCFbjb9Ijim/aozLMRVnUMHn3x+oLMWuZvs3/VfbRESLpOck7Zn231vR\n7sm0z3IaPx5uugnmzoWmpmz0+PTp2R3Hiy9CSwsMGZLNmjtkyLbbO9rX1fZDhmTTu/ef+0Mza0/R\nCaM7/LbSy044AX73uyxpNDVt3f/KK1ni2LRpaxJp3d7Rvueeg6eeyt9+0ybYvDlLGm2Tyd/+Br//\n/dZkUvm9J/t6co5ly7JuvWqybBnMm1d2FNtburT64nJMxSl0pLekE4GmiJiYHl8CRERcUdHmD6nN\nbEkDgbURsXfbtpL+CEyJiNltrlEFHStmZv1PV0d6F32HMRcYJ2kssBaYDJzbps2twPnAbOB/AXek\n/TOA6yV9j6wrahwwp+0FuvqCzcysewpNGKkmcSEwk6zAfk1ELJY0FZgbEb8FrgGuk7QceJYsqRAR\nj0iaDjwCbAYuqPlJo8zMqli/n3zQzMz6Rr8d6S3pGklPS1pYdiytJI2WdIekRZIekvS5KohpZ0mz\nJS1IMU0pO6ZWkgZImi9pRtmxtJL0uKQH089ruy7QMkgaLunGNIB1kaSGkuM5JP185qfvz1XJ//WL\nJD0saaGk6yXtVAUxfT793pX6ftDe+6WkEZJmSloq6U+Shnd2nn6bMIBryQYEVpNXgC9GxBHAG4DP\ntDNQsU9FxMvAWyLiGOBo4AxJ1TKe5fNkXY7VZAvQGBHHRES1/Jx+APw+Ig4HjgIWlxlMRCxLP59j\ngeOAF4BbyoxJ0n7AZ4FjI+L1ZN3tk0uO6QiygcnHk/3unSnpoJLCae/98hLgtog4lKx2fGlnJ+m3\nCSMi/gKsLzuOShHxVEQ8kLb/QfaLXfrYkYjYlDZ3JvtFKr0fUtJo4B3Az8uOpQ1RRb8XkoYBb4qI\nawHSQNbnSw6r0mnAoxFRDVNfDgR2lTQIGAqsKTmew4HZEfFyRLQAfwbOLiOQDt4vKwdN/wI4q7Pz\nVM0vRq2R9Fqyvypm77hl8VLXzwLgKWBWRMwtOybge8C/UgXJq40A/iRprqSPlx0McCDwd0nXpi6g\nn0kaUnZQFc4BflV2EBGxBvgO8ATZIN8NEXFbuVHxMPCm1PUzlOwPpANKjqnS3hHxNGR/7AJ7d3aA\nE0YBJO1GNs3J59OdRqkiYkvqkhoNNKR5ukoj6X8AT6e7MVFdgzVPiojjyX65PyOp7LUNBwHHAj9O\nXUCbyLoSSidpMNkUPjdWQSx7kP3FPBbYD9hN0vvKjCkilgBXALOA3wMLgJYyY+pEp3+8OWH0snQ7\n/Gvguoj477LjqZS6Mu4EJpYcyknAOyU9RvbX6Vsk/Z+SYwIgItam78+Q9cuXXcdYDayKiNax578m\nSyDV4AxgXvpZle004LGIWJe6f24G3lhyTETEtRFxfEQ0AhuAalru7GlJ+wBI2hf4W2cH9PeEUW1/\nnQL8J/BIRPyg7EAAJI1s/fRD6so4HVhSZkwR8eWIGBMRB5EVJu+IiA+WGROApKHp7hBJuwJvI+tW\nKE3qMlgl6ZC061Sq54MC51IF3VHJE8CJknZJs1qfSskfDgCQNCp9HwO8C/hlmeGw7fvlDOBDaft8\noNM/cKtxLqlcJP0SaAT2kvQE2bQh15Yc00nA+4GHUs0ggC9HxB9LDOs1wC/SVPMDgBsi4vclxlPN\n9gFuSdPNDAKuj4iZJccE8DmyWQ8GA48BHy45HlKf/GnAJ8qOBSAi5kj6NVm3z+b0/WflRgXATWky\n1dbBx6V8YKG990vg34EbJX0EWEm2nMSOz+OBe2Zmlkd/75IyM7M+4oRhZma5OGGYmVkuThhmZpaL\nE4aZmeXihGFmZrk4YVi3pCm3P13x+BRJt5YZU3v6Kq40QPI+SfPSeJzK537W2azFkib11czGklak\nsQE9Pc9USW/tpM0pkt7Q02tZdXDCsO4aAVzQZl+1DurpdlxpwGMepwELI+K4iLh7m4tHfCLNK7Qj\nZwFHdCfGtiQN7KRJr/w7RcSUiLijk2aNVMEUHdY7nDCsu74JHJRmUL0i7du9YqGf61obSjpWUnOa\nAfYPrfPXVEqzsf5A0t2S/irp7LR/mzsEST+S9MG0vULSN9IiPnMkHSPpj5KWS6ocgTxc0m8lLZH0\nk4pznS7pHkn3S7ohjV5uPe+/S7ofeE+bOMdKul3ZIkuzlC2adRTZJHOT0s9j5zbH3Cnp2LS9UdLX\nJD2Qrj0q/QX+TuDKdPyBkg5KP6u5ku5qnRok7b83Xf9ySRsrfk5/lvTfwKK075Z0/EOSPlYZUnv/\noCm27ypbhGiWpL3S/qPTNR+QdJO2TjVzbcW/0wpJTekO60FliyyNBT4FfCG9rpMkvSfFs0BSc3tx\nWBWLCH/5q8tfZLOCLqx4fArZfPuvIXtDuofsL8tBwN3AXqnde8nWdm97vmvJpi2BbB2B5RXnnVHR\n7kfAB9P2CuATafu7wANk6yCMBJ6qOH5Tildk68ufDewF3AUMSe3+N/CVivN+qYPXPQM4L21/GLgl\nbZ8P/LCDY+4kW9gHsgWa3pG2ryCbOqb19Z9dccxtwOvS9gTg9rR9K/DetP1J4PmK17kRGFNxjj3S\n912Ah4ARFa9vz3bi3AJMTttfbX09wIPAyWl7KvDdtjGnc16Qtj8N/CxtTyFbVKz1GguB16TtYWX/\nP/ZX17767VxSVpXmRJrtVdIDwGuB54AjgVmSWhcn6mhhm98ARMRiSZ3OzZ+03n08BOwa2WJRmyS9\npGzxoda4Vqa4fgWcDLwMjAfuTnENJktyrW7o4HpvIJtEDuA6sjf9rng5ts7lNY+sK2sbyiY+fCPZ\nPD+tdwODK64/KW3/EvhWxaFzIuKJisdfkNS6KM5o4GBgR8vOtgDT0/b/JZsHaRgwPLIFeCBbaGd6\newezddW9eWz9GbX1F7K5zaaTzShr/YgThvWmlyu2W8j+fwl4OCJOav+QDo9vfaN8hW27Tnfp4Jgt\nbY7fwtb/32377COdf2ZEvL+DWF7oYH9P+/83V2y3/ozaGgCsj2z9ix1dv23X0qsxSzoFeCvQEBEv\nS7qT7X92nWm9Vt4ZoVt//h29LiLiAkknAGcC8yQdGxFVtXKmdcw1DOuujcDuOdotBUZJOhGy9UKU\nbwGn1jeplcB4SYOVLZJzas74Kt/kGlLtYQDZCnF/Ae4DTpL0uhTXUEkH5zjvPWTTegOcB/y/nPG0\nF1eljcAwgIjYCKyQ9Gr9RNLr0+Z9bK2r7GjN6uFkSedlZZ++OjFHbAMrzv1+4C+Rza66Tls/+fUB\nsq68vF59XZDVYCJibkRMIVt/oZpWoLNOOGFYt0TEOrLunIXaWvTepklqt5nsTeiK1E21gKxbpd32\n7Ry/mqwL5GFgGjB/B8d0dL45wFVkxeBHI+KWiPg72VoAv5L0IFkiODTHeT8HfDi9lvcDn99B2/Zi\n6ejc04B/TUXjA9O5P5oKzQ+TFcUBLgK+mK7/OrIuv/b8ERgsaRHwDeDeHDG8AEyQ9BDZp5v+Le0/\nH/h2uuZRFfvzvK5bgXe1Fr2Bb6X/MwuBuyNiYQfHWRXy9OZm/YikIRHxYto+h6xI3VG9oKvn3hgR\nee4arU65hmHWvxwn6Sqyrq31wEd68dz+69F2yHcYZmaWi2sYZmaWixOGmZnl4oRhZma5OGGYmVku\nThhmZpaLE4aZmeXy/wEK4DAec3DG+wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efbe4b44fd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "class Example2():\n",
    "    def __init__(self, h=0.05):\n",
    "        self.h = h\n",
    "        self.q = 0\n",
    "        self.points = np.array([\n",
    "                [1.0, 1.0],\n",
    "                [2.0, 1.0],\n",
    "                [2.0, 2.0],\n",
    "                [1.0, 2.0],\n",
    "                [0.5, 1.5]\n",
    "            ])\n",
    "    def __call__(self, p):\n",
    "        x = p[:, 0]\n",
    "        y = p[:, 1]\n",
    "        return np.exp(x/y)\n",
    "    \n",
    "f = Example2()\n",
    "fe = exact_quad(f)\n",
    "print(fe)\n",
    "maxit = 10\n",
    "error = np.zeros(maxit, dtype=np.float)\n",
    "for n in range(1, maxit+1):\n",
    "    fn = polygon_quad(f, n)\n",
    "    error[n-1] = np.abs(fn-fe)\n",
    "print(error)\n",
    "print(error[:-1]/error[1:])\n",
    "\n",
    "g = plt.figure()\n",
    "axes = g.gca()\n",
    "axes.set_xlim([0.5,2.1])  \n",
    "axes.set_ylim([0.9,2.1])\n",
    "axes.set_axis_off()\n",
    "coord = [[1,1], [2,1], [2,2], [1,2], [0.5,1.5]]\n",
    "coord.append(coord[0])\n",
    "xs, ys = zip(*coord)\n",
    "plt.plot(xs, ys,linewidth=1.0)\n",
    "plt.title('the figure of Integral area')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "x = np.arange(1,maxit+1)\n",
    "y = error\n",
    "plt.plot(x,y)\n",
    "xlabel('the number of integral points')\n",
    "ylabel('error')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "78.2083333333\n",
      "[  2.08333333e-01   5.68434189e-13   5.68434189e-13]\n",
      "[  3.66503876e+11   1.00000000e+00]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWcAAAEKCAYAAADO0pQJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFndJREFUeJzt3XuwXWV5x/HvkxATjIQMY6oJJRBqgDYFM0yLqaQalIvU\nKkixWEFULp1g44Qi5HipFdGqgWKhhXYiZagyAxQZMENFgXgh4g0sGS7eUi4G1MglUUDQaMjbP9Y6\nnp3DuZ+997vWXt/PzB7Ovpy1nrM3+e3nrPXs90RKCUlStUzJXYAk6fkMZ0mqIMNZkirIcJakCjKc\nJamCDGdJqiDDueYiYu+I2BERbXktI+KMiPhZRDwVEXtExNMRsU87tt1JEbFfRGyIiCcjYkXueiYj\nIt4eEV/LXYfyMpxrJiIeiojXDLq5LcPqEbELcCFweEppVkppa0ppt5TSj9qx/Q5bBXw5pbR7SumS\nwXdGxFci4pSxbCgiroiI89pe4fj4AYSGM5zV6qXAdOD7nd5RRExt8yb3Br7b5m22XQd+7tH257/x\nmvKFq5GI+AwwH7ixPOxwdv9dwEkRsSkiHouI97d8T0TEeyPi/oh4PCKuiYjZQ2x7IfCD8urPI2Jd\nefuOiNi3/HqPiLixPHTw7Yj4SP+v30MdXmntVstf1W+PiE9GxBPAh8rbT4mI70XEloj4QkTMH+Hn\nf2NE3BcRWyPiyxGxf3n7l4DDgEvL5+VlozyPr46IRyLirIh4NCJ+EhHvKO87HTgRWFVua215+9yI\nuK58fh+IiHe3bG9GRHy6rOu7EXFORDzScv9DEbEqIu4GfhkRUyKir3xNnip/pmNHqnlQ/ddGxOaI\n+HlEfDUi/qjlvisi4t8j4vMR8TSwLCJeEBH/XP7/sbm8f3r5+Nnla/pY+RrcGBHzxlqLOiil5KVG\nF+Ah4LCW63sDO4A1wAuAg4BfA/uX968EvgHMBaYB/wFcNcy29waeA6LltueAfcuvrwGuouiu/xB4\nGFg/6HuntHzvV4BTyq/fDvwWeBdFUzAdOAbYCOxX3vZ+4OvD1LYf8EvgNcBU4Bzg/4BdBu9rmO9v\nreXVZS0fKrd1NPAMsHt5/xXAeS3fG8B3gA+Uj98HuB84orz/E+X2ZwHzgLuBhwe9ZneV900vb/sr\n4CXl128uf7aXtDxX60f4Wd4BvLB8PT8JbGi57wrg58CS8vp04F+AzwG7AzOBtcA/lffvAbypfNxM\n4L+B63P/f+4lGc51u5T/0F/Tcr0/FOe23PZt4K/Lr7/HzmE+F/hNa4gOsa3WgN0B7FuG52+Al7Xc\n9xHGF84/GrS/m4B3tlyfUobkXkPU9g/ANS3XA/gx8KrB+xrmeRsczs8MqvVR4JDy68HhfMgQtb8X\nuLz8+gGK4/T99506RDi/fZTXdQPwhpbnathwHvR9s8vXaLeW2v9r0GN+CSxouf5nwIPDbG8xsCX3\n/+deErugXvFoy9fPAi8qv94buCEidpTXg6JrfAmweRzbn0PRNf645bZHhnnscAY/fm/g4oi4sKW2\nBOw5xGPnAZv6r6SUUnnoYM9x1tBvS0ppR8v11udssL2BPSNia0udU4D1LbWN9ry03k9EnAz8PUUX\nDkXX+uLRii4PG30MOL58fCovLwaeHrz/iJhD0WX/b0T03zyl/BmIiF2Bi4CjKII+gBdFRKQyrZWH\n4Vw/4/0H8zBFx/jNSe73cWA78PsUv9ID7NVy/zPlf19I0alBcYKx1eDaHwY+mlK6egz7/ynwx4Nu\n24tBodcmg+t8hKLT3H+Yx/+U4nnpP2Y/1HHz322zPK7+KYrfaL5Z3raBMjBH8VbgDRS/PT0cEbtT\nHMZo/d7W+p+geONZlFIa6s34PcBC4E9TSo9HxMspDsH0v1EqE08I1s/PKA4ztBrpH/Ua4GP9J9oi\nYk5EvHGExw+5rbLLvB44NyJ2jYgDgJNb7n8C+AnFickp5YnAPxjlZ1kDvL//hFZE7B4Rxw/z2GuB\n10fEYRGxS3ky9NfAZN90hvIoOz/HdwBPlyf1ZkTE1IhYFBF/Ut7/WeB95cm1PYG/G2X7MykORTxR\nPlfv5PlvPMPZDdhGcdJ2JvBxRgjRsvu9DLio7KKJiD0j4siW7f0KeCoi9gDOHWMd6jDDuX4+AXyw\nnAw4q7xt8D/O1usXU5wAuiUinqQ4OXjICNsfaVvvpvjVdzPwaYqTg9ta7j+dYt74CYoThl8f6QdJ\nKX2u/HmuiYhfAPcArxvmsRuBk4BLKLr411Mco90+TN2j/Vwj3X85sKh8jq8v35j+kuJ47EPAYxSB\nN6t8/HkUb0wPAbdQhHXr87LTvlNK36eYJ/8WxZvtIuD2Uerr9xmK3zh+AtxH8XqOpo/it51vlc/z\nLRQnWKE4pPFCitfsGxTnAVQB4WElTVREfIJiwuCduWupkohYDpyQUjosdy2qLztnjVlE7B8RB5Zf\nH0IxlXB93qryi4iXRsQro7A/xXHcxj8vmhxPCGo8dgOujoi5FMdlL0gp3Zi5pip4AcXx832AXwBX\nU8yTSxPmYQ1JqiAPa0hSBRnOklRBhrMkVZDhLEkV5LSGpOeJYBbFMqyHAv+T0u/WEVGX2DlLIoJp\nERwawbkRfJ3iE4grKFYi/HQEN0aM+SPmagNH6aQGiiCA/YEjysurgAeBdcCtwO0p8avysdOB5RTr\nbX8B+MeUeDhH3U1iOEsNEcEc4HCKMD68vPnW8vKllHh8lO+fBZxNsbDTFcDHU2JL5ypuNsNZ6lER\n7AosZaA7XgDcxkAgb0xp/MuCRjAX+CDFX3C5EPjXlHi2XXWrYDhLPSKCKRQr5/WH8Sso/mRW/6GK\nO1Lit23c337AR4FXAh8GrkiJ7SN/l8bKcJZqLIL5DITxa4EtDHTGX02Jp7pQwyEUS7/OA94HfG4i\nHbl2ZjhLNRLB7sAyBgJ5DwY643W5TtSVJxiPBFZTLN7f5/jd5BjOUoVFMI3i8ET/SbyDKBbp7++O\n706JHcNvobvKQyt/Q3G44z7gfSlxX96q6slwlipkPCNuVeb43eQZzlJmkx1xqzLH7ybOcJa6rFMj\nblXm+N34Gc5Sh3V7xK3KHL8bO8NZ6oAqjLhVWTl+txqYi+N3QzKcpTao6ohblZUnP4+imJF2/G4Q\nw1magLqNuFWZ43dDM5ylMeiVEbcqc/xuZ4azNIwIfo/ieHF/ICd6ZMStyhy/KxjOUqmJI25VNmj8\n7pPAxU0avzOc1VjDjLjdw0AYN2bErcqaOn5nOKtRRhhxW0cx4vZkxvI0gqaN3xnO6mmOuPWWJo3f\nGc7qKY64NUMTxu8MZ9WaI27N1svjd4azascRNw1WHr46G3gXPTJ+Zzir8soRtz9n4FCFI24aUi+N\n3xnOqhxH3DRZvTB+ZzirEhxxUyfUefzOcFYWjripW+o6fmc4qysccVNudRu/M5zVEUOMuL2aYsSt\nP4wdcVMWdRm/M5zVNiOMuK2jGHF7LGN50k6qPn5nOGvCHHFTL6jq+J3hrDFzxE29rGrjd4azRuSI\nm5qmKuN3hrN2Uh6HO4ziMIUjbmqklvG71cCzZBi/M5wbbtCI2xHAgTjiJgG/O5T3VuAjwHcpxu/u\n7cq+DedmccRNGr8c43eGcwM44ia1RzfH7wznHuSIm9RZ3Ri/M5x7gCNuUh6dHL8znGtqiBG3rQyE\nsSNuUhd1YvzOcK4JR9ykamv3+J3hXHERLACuBF6OI25S5Q0av/tKSpwyoe0YztUWwX8CTwEfcMRN\nqo8I9qI43zN3It+/S5vrURtFMA84DlhoMEu1M6kTg1PaVYU64kzgyiotYyipO+ycKyqC2cCpwMG5\na5HUfXbO1bUcuCklNuUuRFL32TlXUAQzgJUUYzmSGsjOuZreBmxIiXtyFyIpDzvniolgKnAOcHru\nWiTlY+dcPcdSfBS7qwt7S6oWw7lCyo9/9gGrXTVOajbDuVqWAbOAtZnrkJSZ4VwtfcAFrpkhyROC\nFRHBYuAg4JjctUjKz865Os4BLkqJbbkLkZSf4VwBEewDvA5Yk7kUSRVhOFfDe4DL/Oslkvp5zDmz\nCOYAJwKLctciqTrsnPNbAVyXEptzFyKpOuycM4pgJnAGsDR3LZKqxc45r9OA9SmxMXchkqrFzjmT\nCKYBZwFvzl2LpOqxc87nBODBlLgjdyGSqsfOOYNygaNV5UWSnsfOOY+jgR3AzbkLkVRNhnMefcD5\nLgsqaTiGc5dFsASYD1ybuxZJ1WU4d18fcGFKbM9diKTq8oRgF0VwAHAocFLuWiRVm51zd50NXJoS\nz+QuRFK12Tl3SQTzgOOAhblrkVR9ds7dcyZwZUpsyV2IpOqzc+6CCGYDpwIH565FUj3YOXfHcuCm\nlNiUuxBJ9WDn3GERzABWAkfmrkVSfdg5d97bgA0pcW/uQiTVh51zB0UwleKvap+euxZJ9WLn3FnH\nAluB9bkLkVQvhnOHlMuC9gGrXeBI0ngZzp2zDJgFrM1ch6QaMpw7pw+4ICV25C5EUv14QrADIlgM\nHAgck7sWSfVk59wZ5wAXp8S23IVIqifDuc0i2Ad4HbAmcymSasxwbr/3AJelxJO5C5FUXx5zbqMI\n5gAnAoty1yKp3uyc22sFcF1KbM5diKR6s3NukwhmAmcAS3PXIqn+7Jzb5zRgfUpszF2IpPqzc26D\nCKYBZwFvzl2LpN5g59weJwAPpsQduQuR1BvsnCepXOBoVXmRpLawc568o4EdwM25C5HUOwznyesD\nzndZUEntZDhPQgRLgPnAtblrkdRbDOfJ6QMuTIntuQuR1Fs8IThBERwAHAqclLsWSb3HznnizgYu\nTYlnchciqffYOU9ABPOA44CFuWuR1JvsnCfmTODKlNiSuxBJvcnOeZwimA2cChycuxZJvcvOefyW\nAzelxKbchUjqXXbO4xDBDGAlcGTuWiT1Njvn8TkZ2JAS9+YuRFJvs3MeowimUozPnZ67Fkm9z855\n7I4FtgLrcxciqfcZzmNQLgvaB6x2gSNJ3WA4j80yYBawNnMdkhrCcB6bPuCClNiRuxBJzeAJwVFE\nsBg4EDgmdy2SmsPOeXSrgItSYlvuQiQ1h+E8gggWAEcBn8pdi6RmMZxHdhZwWUo8mbsQSc3iMedh\nRDAHOBFYlLsWSc1j5zy8FcB1KbE5dyGSmsfOeQgRzATOAJbmrkVSM9k5D+00YH1KbMxdiKRmsnMe\nJIJpFCcCj89di6TmsnN+vhOAB1PiztyFSGouO+cW5QJHq8qLJGVj57yzo4EdwM25C5HUbIbzzvqA\n810WVFJuhnMpgiXAfODa3LVIkuE8oA+4MCW25y5EkjwhCERwAHAoxce1JSk7O+fC2cClKfFs7kIk\nCeyciWAecBywMHctktTPzhnOBK5MiS25C5Gkfo3unCOYDZwKHJy7Fklq1fTOeTlwU0psyl2IJLVq\nbOccwQxgJXBk7lokabAmd84nA3elxL25C5GkwRrZOUcwlWJ87vTctUjSUJraOR8LbAXW5y5EkobS\nuHAulwXtA1a7wJGkqmpcOAPLgFnA2sx1SNKwmhjOfcAFKbEjdyGSNJxGnRCMYDFwIHBM7lokaSRN\n65xXARelxLbchUjSSBoTzhEsAI4CPpW7FkkaTWPCGTgLuCwlnsxdiCSNphHHnCOYQ7GQ/qLctUjS\nWDSlc14BXJcSm3MXIklj0fOdcwQzgTOApblrkaSxakLnfBqwPiU25i5EksaqpzvnCKZRnAg8Pnct\nkjQevd45nwA8mBJ35i5EksajZzvncoGjVeVFkmqllzvno4EdwM25C5Gk8erlcO4DzndZUEl11JPh\nHMESYD5wbe5aJGkiejKcKbrmC1Nie+5CJGkieu6EYAQHAIdSfFxbkmqpFzvnc4BLU+LZ3IVI0kT1\nVOccwTzgTcDC3LVI0mT0Wud8JnBlSmzJXYgkTUbPdM4RzAZOBQ7OXYskTVYvdc7LgZtSYlPuQiRp\nsnqic45gBrASODJ3LZLUDr3SOZ8M3JUS9+YuRJLaofadcwRTKcbnTstdiyS1Sy90zscCW4D1uQuR\npHapdTiXy4L2Aatd4EhSL6l1OAPLgFnA2sx1SFJb1T2c+4ALUmJH7kIkqZ1qe0IwgsXAgcAxuWuR\npHarc+e8CrgoJbblLkSS2q2W4RzBAuAoYE3uWiSpE2oZzsBZwGUp8VTuQiSpE2p3zDmCORQL6S/K\nXYsktSrHe/cDjgD+AvjVRLdVu3AGVgDXpcTm3IVIUtkwHk4RyIeXN98KfAZYN+HtplSfz25EMBN4\nCFiaEhtz1yOpeSLYFVhKEcZHAAuA2ygC+VZgYzs+FFe3zvk0YL3BLKlbIpgCvJyBMF4C3E3RFa8A\n7kiJ37Z9v3XpnCOYBtwPHJ8Sd+auR1LvimA+A2H8Wor1e/o74692YxihTp3zW4AHDGZJ7RbBLOAw\nBgJ5D4rO+BZgVUo83PWa6tA5l2dA7wHOSYkv5q5HUr2Vv4kfwkAYHwR8i4Hu+O7cy0LUpXM+GngO\nuDl3IZLqZ9CI2xHAq4EHKbrjc4HbU5r42Fsn1KVzvg1YkxJX5a5FUj2MMOJ2K/CllHg8V21jUflw\njmAJcDWwMCW2565HUjV1a8StW+oQzjdQvMtdkrsWSdUxzIjbPQyEcUdG3Lql0uEcwQEU73wLUuLZ\n3PVIymuEEbd1FCNuT2Ysr62qHs6XA5tS4rzctUjqvhFG3G4F1uUYceuWyoZzBPOA+yiONW/JXY+k\nzqvDiFu3VDmczwemp8TK3LVI6oxyxG1/BsL4VQyMuN1KBUfcuqWS4RzBbOAB4OCU2JS7HkntU/cR\nt26p6odQlgM3GcxS/Y0y4raamo24dUvlOucIZlAsC3pkStybux5J49PrI27dUsXO+WTgLoNZqo8R\nRtwuoVhJsmdG3LqlUp1zBFOBHwCnpsT63PVIGlqTR9y6pWqd87EU77hfy12IpAGjjLi9hQaNuHVL\nZTrncqTm28DHU+KG3PVITTbEiFv/Km79x40bO+LWLVXqnJcBs4C1meuQGmmEEbergdNS4rFctTVR\nlTrnLwKfTYnLc9ciNUGvreLWayoRzhEsBj4P7JsS23LXI/UiR9zqpSrhfBWwISUuyF2L1EuGGHHb\nys5/qNQRt4rKHs4RLAC+Q7EsaMf/oq3Uyxxx6x1VCOd/A55JifdmLUSqIVdx611Zw7k8O/xDYFFK\nbM5WiFQTjrg1R+5w/jAwNyX+NlsRUsWNMOK2jmIVN0fcelC2cI5gJsUCR0tTYmOWIqQKcsRNkPdD\nKKcB6w1mNd0oI24rcMStkbJ0zuVJjPspVqu6s+sFSJk54qbR5Oqc3wI8YDCrKUYYcbsFWOWImwbr\neudcnm2+BzgnJb7Y1Z1LXeKImyYrR+d8NPAccHOGfUsdMcqI27k44qZxytE53wasSYmrurpjqc2G\nGHELdv5DpY64acK6Gs4RLKFYfnBhSmzv2o6lNhhlxG0d8ENH3NQu3Q7nGyg6iku6tlNpgsoRt8UM\nhPErcBU3dUnXwjmCAyi6jAUp8WxXdiqNkyNuqopuhvPlwKaUOK8rO5TGwFXcVFVdCecI5gH3URxr\n3tLxHUrDcMRNddGtcD4fmJ4SKzu+M6mFq7iprjoezmWn8pvyqmey1W0B/BhH3FQz3eqcXwSeBFQ2\nyRE31U32v4QiSXq+KbkLkCQ9n+EsSRVkOEtSBRnOklRBhrMkVZDhLEkVZDhLUgUZzpJUQYazJFWQ\n4SxJFWQ4S1IFGc6SVEGGsyRVkOEsSRVkOEtSBRnOklRBhrMkVZDhLEkVZDhLUgUZzpJUQYazJFWQ\n4SxJFWQ4S1IF/T+64bs2r2C0jgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efc0e61e7f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEPCAYAAABV6CMBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGtxJREFUeJzt3XnQHPV95/H3R8IQMCCDwWuKy8byBakYhBEkAusxGCO7\nHIS1xBaFAvKypYCLSuJgL0k2i6TElxzjxBhTXmKWMngdYYdL4EDEYj2wiEMHt9CFrRUCZFMxAstc\nUqTv/tG/kVrDzPNMj6bn/LyqptTT093z62bQV/379adbEYGZmVmjxnS6AWZm1ltcOMzMrBAXDjMz\nK8SFw8zMCnHhMDOzQlw4zMyskNILh6QpklZJWiPp0hqff0HSCkmPSrpL0uG5z7ZJeljSI5JuKbut\nZmY2OpWZ45A0BlgDnAY8DywFpkfEqtwyk4GHIuJ1SRcCQxExPX32m4jYv7QGmplZYWWfcUwE1kbE\n+ojYCswHpuYXiIh7IuL19PZB4NDcxyq5fWZmVlDZheNQYEPu/bPsWhiqXQDckXu/l6Qlku6XNLXe\nSmZm1j57dLoBFZJmAMcDk3Ozj4yIjZLeDfxM0uMRsa4zLTQzMyi/cDwHHJF7f1iatwtJHwP+CvhI\n6tICICI2pj/XSRoGjgPWVa3rm22ZmTUhIpoaDii7q2opMF7SkZL2BKYDC/ILSDoO+B5wZkT8Ojf/\nbWkdJB0E/AHwVK0viQi/WvSaPXt2x9vQTy8fTx/Pbn3tjlLPOCJim6SLgYVkReqaiFgpaS6wNCJu\nB74BvBX4iSQB6yPiLOCDwP+UtC2t+7XIXY1lZmadUfoYR0TcCby/at7s3PTpddZ7APi9cltnZmZF\nOTluuxgaGup0E/qKj2dr+Xh2h1IDgO0gKXp9H8zM2k0S0aWD42Zm1mdcOMzMrBAXDjMzK8SFw8zM\nCnHhMDOzQlw4zMysEBcOMzMrxIXDzMwKceEwM7NCXDjMzKwQFw4zMyvEhcPMzApx4TAzs0JcOMzM\nrBAXDjMzK8SFw8zMCnHhMDOzQlw4zMysEBcOMzMrxIXDzMwKceEwM7NCXDjMzKwQFw4zMyvEhcPM\nzApx4TAzs0JcOMzMrBAXDjMzK8SFw8zMCnHhMDOzQlw4zMyskL4oHJs3d7oFZmaDoy8Kx7RpsGVL\np1thZjYY+qJw7LsvzJwJ27d3uiVmZv2vLwrHj34EGzbAl77U6ZaYmfW/0guHpCmSVklaI+nSGp9/\nQdIKSY9KukvS4bnPzk/rrZZ0Xr3v2HtvWLAA/u3f4JvfLGtPzMwMQBFR3salMcAa4DTgeWApMD0i\nVuWWmQw8FBGvS7oQGIqI6ZIOAJYBEwABy4EJEfFy1XdEZR82bICTT4avfAVmzChtt8zMep4kIkLN\nrFv2GcdEYG1ErI+IrcB8YGp+gYi4JyJeT28fBA5N02cACyPi5Yh4CVgITBnpyw4/HO64Ay65JDv7\nMDOz1iu7cBwKbMi9f5adhaGWC4A76qz73CjrAnD00XDTTdkZx7JlBVtrZmaj2qPTDaiQNAM4Hphc\ndN05c+bsmB4aGmJoaIjvfx/OPBPuvRfGj29dO83MetHw8DDDw8Mt2VbZYxwnAXMiYkp6/5dARMS8\nquU+Bnwb+EhE/DrNm0423nFhev89YFFE3FC1btTbh6uvhnnzYPFieOc7W7xzZmY9bHfGOMouHGOB\n1WSD4xuBJcA5EbEyt8xxwE+AMyLi57n5+cHxMWn6+DTekf+OuoUDYO5cuPVWuOce2G+/lu2amVlP\n69rCAdnluGRnE2OAayLi65LmAksj4nZJdwG/S1ZYBKyPiLPSujOB/w4E8OWIuK7G9kcsHBFw0UXw\n85/DT38Ke+7Z4h00M+tBXV04yjZa4QDYtg3OPjvLe/zwhzCmL2KPZmbN6+bLcbvC2LFOl5uZtcpA\nFA5wutzMrFW65nLcdjjggCwgePLJ2VVWTpebmRU3UIUDdqbLP/pROPhgOOOMTrfIzKy3DExXVZ7T\n5WZmzRvIwgEwaRI70uVPP93p1piZ9Y6B66rKmzoVfvWrrLvK6XIzs8YMdOEAmDULNm6ET37S6XIz\ns0YMRABwNE6Xm9mgcXK8BfvgdLmZDRInx1vA6XIzs8a4cOQ4XW5mNrqBHxyv5nS5mdnIXDhqcLrc\nzKw+d1XV4XS5mVltLhwjcLrczOzN3FU1CqfLzcx25cLRAKfLzcx2cgCwQU6Xm1k/cXK8TfvgdLmZ\n9Qsnx9vE6XIzMxeOwpwuN7NB58HxJjhdbmaDzIWjSU6Xm9mgclfVbnC63MwGkQvHbnK63MwGjbuq\nWsDpcjMbJC4cLeJ0uZkNCgcAW8jpcjPrFU6Od9E+OF1uZr3AyfEu4nS5mfU7F44SOF1uZv3Mg+Ml\ncbrczPqVC0eJnC43s37krqqSOV1uZv2m9MIhaYqkVZLWSLq0xuenSFouaaukaVWfbZP0sKRHJN1S\ndlvL4nS5mfWTUruqJI0BrgROA54Hlkq6NSJW5RZbD5wPfLHGJl6JiAlltrFdnC43s35R9hjHRGBt\nRKwHkDQfmArsKBwR8Uz6rFYYo6lrjLuV0+Vm1g/K7qo6FNiQe/9smteovSQtkXS/pKmtbVpnXHYZ\nTJwI06bBli2dbo2ZWXHdPjh+ZERMBM4F/lHSuzvdoN0lwXe/C/vuCzNnwvbtnW6RmVkxZXdVPQcc\nkXt/WJrXkIjYmP5cJ2kYOA5YV73cnDlzdkwPDQ0xNDTUVGPbpZIu//jHs3T55Zd3ukVm1u+Gh4cZ\nHh5uybZKvVeVpLHAarLB8Y3AEuCciFhZY9lrgdsj4sb0/m3AqxGxRdJBwGJgatXAetfdq6qITZvg\nlFOyM48v1ro0wMysJF17r6qI2AZcDCwEVgDzI2KlpLmSPgUg6cOSNgBnA9+T9ERa/YPAMkmPAHcD\nX6suGr2uki7/zneyGyKamfUC3x23Czz1VJYuv+46p8vNrD269ozDGuN0uZn1EheOLuF0uZn1Ct/k\nsIs4XW5mvcCFo8s4XW5m3c6D413Izy43s7L5meM9vg+1+NnlZlYmX1XVh/zscjPrVi4cXczPLjez\nbjRq4ZA0VtIX2tEYe7NKuvyKK5wuN7PuMGrhSLcNOacNbbE6Ks8uv+SS7OzDzKyTGhocl/QPwFuA\nG4BXKvMj4uHymtaYfh0cr2XxYjjrrKyIfPjDnW6NmfWy0q+qkrSoxuyIiFOb+dJWGqTCAXDrrdml\nuvfeC+PHd7o1ZtardqdwNBQAjIiPNrNxaz2ny82s0xq6qkrSOEnfkrQsvS6XNK7sxllts2bBeedl\n6fLNmzvdGjMbNI12Vd0IPAn8IM36Y+BDETGtxLY1ZNC6qioi4MIL4Re/cLrczIprxxjHoxFx7Gjz\nOmFQCwc4XW5mzWtHcvw1SSfnvnAS8FozX2it43S5mXVCo3fHvRC4LjeusQk4v5wmWRGVdPkpp8Ah\nh/jZ5WZWvlELh6QxwPsj4kOS9geIiN+U3jJrWCVdPmlSdpXVjBmdbpGZ9bNGxziWRURXRs4GeYyj\n2ooVcOqpfna5mY2uHYPjXwf+nTcnx19s5ktbyYVjV06Xm1kj2lE41tWYHRFxVDNf2kouHG/mdLmZ\njabU5Hga45gREYub+QJrP6fLzaxMjdwddztwZRvaYi3kdLmZlaXRHMfdkv6zpKZOa6wzLrsMTjgB\npk2DLVs63Roz6xeNjnFsBvYBtgGvAyIb49i/3OaNzmMcI3O63MxqaUdyfBwwE/hyKhbHAKc384XW\nXk6Xm1mrNVo4vgucxM4nAW7G4x49w88uN7NWavSWIydGxARJjwBExCZJvh9rD3G63MxapdHCsVXS\nWCAAJB0MbC+tVVaKyrPLTz0VDj7Y6XIza06jXVVXADcD75D0FeA+4KultcpKc8wxcNNN2RnHsmWd\nbo2Z9aKGrqoCkPQB4DSyK6rujoiVZTasUb6qqjlOl5sNttKfOQ4QEauAVc18iXUfp8vNrFkNFw7r\nP7NmwcaNWbr8nntgv/063SIz6wUNd1V1K3dV7R4/u9xsMJV+d9xu5sKx+5wuNxs87UiON03SFEmr\nJK2RdGmNz0+RtFzSVknTqj47P623WtJ5Zbd1UDldbmZFlHrGkW7JvobsaqzngaXA9DTQXlnmCGB/\n4IvAgoi4Kc0/AFgGTCC7kms5MCEiXq76Dp9xtMimTdmzy2fO9LPLzfpdN59xTATWRsT6iNgKzAem\n5heIiGci4klSuDDnDGBhRLwcES8BC4EpJbd3oFXS5VdckXVZmZnVUvZVVYcCG3LvnyUrJs2s+1ya\nZyVyutzMRuNhUHsTp8vNbCRln3E8BxyRe39YmtfoukNV6y6qteCcOXN2TA8NDTE0NFRrMStg0iT4\n/vfhzDOdLjfrB8PDwwwPD7dkW2UPjo8FVpMNjm8ElgDn1LpdiaRrgdsj4sb0Pj84PiZNH5/GO/Lr\neXC8RFdfDfPmOV1u1m+6dnA8IrYBF5MNbK8A5kfESklzJX0KQNKHJW0Azga+J+mJtO4m4O/ICsZD\nwNzqomHl87PLzayaA4A2KqfLzfqPk+M9vg+9wOlys/7StV1V1j+cLjezChcOa5ifXW5m4NuqW0F+\ndrmZuXBYYU6Xmw02d1VZU5wuNxtcLhzWtHy6/OmnO90aM2sXd1XZbvGzy80GjwuH7TY/u9xssDgA\naC3hdLlZb3FyvMf3oV84XW7WO5wct67gdLnZYHDhsJZyutys/3lw3FrO6XKz/ubCYaVwutysf7mr\nykrjdLlZf3LhsFI5XW7Wf9xVZaVzutysv7hwWFs4XW7WPxwAtLZxutysezg53uP7MEicLjfrDk6O\nW89wutys97lwWNs5XW7W2zw4bh3hdLlZ73LhsI5xutysN7mryjrK6XKz3uPCYR03aRL80z/BH/6h\n0+VmvcBdVdYVzjoLXnjB6XKzXuDCYV3D6XKz3uAAoHUVp8vN2sPJ8R7fB9uV0+Vm5XNy3PpKJV3+\nzDNOl5t1IxcO60qVdPmddzpdbtZtPDhuXevAA7PC4XS5WXdx4bCu5nS5WfdxV5V1PafLzbpL6YVD\n0hRJqyStkXRpjc/3lDRf0lpJD0g6Is0/UtKrkh5Or6vKbqt1L6fLzbpHqV1VksYAVwKnAc8DSyXd\nGhGrcotdALwYEe+V9FngG8D09NnTETGhzDZa73C63Kw7lH3GMRFYGxHrI2IrMB+YWrXMVOAHafpf\nyIpMRVPXGFv/mjULzjsvS5dv3tzp1pgNprILx6HAhtz7Z9O8mstExDbgJUkHps/eJWm5pEWSTi65\nrdYjLrsMTjgBpk2DLVs63RqzwdONV1VVzjI2AkdExCZJE4BbJB0dEb+tXmHOnDk7poeGhhgaGmpH\nO61DJLjqqixdPnOm0+VmjRgeHmZ4eLgl2yr1liOSTgLmRMSU9P4vgYiIebll7kjLPCRpLLAxIt5R\nY1uLgEsi4uGq+b7lyIB67TU4/XQ48US4/PJOt8ast3TzLUeWAuPTFVJ7kg16L6ha5jbg/DT9R8DP\nACQdlAbXkXQUMB74RcnttR7idLlZZ5TaVRUR2yRdDCwkK1LXRMRKSXOBpRFxO3ANcL2ktcCv2XlF\n1UeAv5W0BdgO/ElEvFRme633OF1u1n6+O671hRUrsnT5ddc5XW7WiG7uqjJrC6fLzdrHhcP6htPl\nZu3RjZfjmjXN6XKz8rlwWN/xs8vNyuXBcetLfna52cj8zPEe3wcrh59dblafr6oyq8HPLjcrhwuH\n9TWny81az4Pj1vecLjdrLRcOGwh+drlZ67irygaG0+VmreHCYQPF6XKz3eeuKhs4Tpeb7R4XDhtI\nTpebNc8BQBtYTpfbIHNyvMf3wTrH6XIbVE6OmzXJ6XKz4lw4bOA5XW5WjAfHzXC63KwIFw6zxOly\ns8a4q8osx+lys9G5cJhVcbrcbGTuqjKrwelys/pcOMzqcLrcrDYHAM1G4HS59Ssnx3t8H6y7OV1u\n/cjJcbMSOV1utisXDrMGOF1utpMHx80a5HS5WcaFw6wAp8vN3FVlVpjT5TboXDjMmuB0uQ0yd1WZ\nNcnpchtULhxmu8HpchtEDgCa7Sany60XdXUAUNIUSaskrZF0aY3P95Q0X9JaSQ9IOiL32V+l+Ssl\nfbzstpo1Q4KrroJ994WZM2H79k63yKxcpRYOSWOAK4EzgGOAcyR9oGqxC4AXI+K9wD8C30jrHg18\nBvgg8AngKklNVUdr3PDwcKeb0JPqpct9PFvLx7M7lH3GMRFYGxHrI2IrMB+YWrXMVOAHafpfgFPT\n9JnA/Ij4j4j4f8DatD0rkf/HbF6tdLmPZ2v5eHaHsgfHDwU25N4/y5v/8t+xTERsk/SypAPT/Ady\nyz2X5pl1rep0uVk/6sarqtwdZT0tny4fNw6WL+90i/rH6tU+nt2g7MLxHHBE7v1haV7es8DhwPOS\nxgL7R8SLkp5L80daF8iuDrDWmTt3bqeb0DdeeAHWrvXxbCUfz84ru3AsBcZLOhLYCEwHzqla5jbg\nfOAh4I+An6X5C4D/LekfyLqoxgNLqr+g2cvJzMysOaUWjjRmcTGwkGwg/pqIWClpLrA0Im4HrgGu\nl7QW+DVZcSEinpL0Y+ApYCvweQc2zMw6r+cDgGZm1l49c5NDSddI+pWkx0dY5ooUGHxU0rHtbF8v\nGe1YSpos6SVJD6fX37S7jb1E0mGSfiZphaQnJP1pneX8+xxFI8fSv8/GSdpL0kOSHknHc3aNZeqG\nsOuKiJ54AScDxwKP1/n8E8BP0/SJwIOdbnO3vho4lpOBBZ1uZ6+8gHcCx6bpfYHVwAeqlvHvs3XH\n0r/PYsd0n/TnWOBBYGLV5xcBV6Xpz5Ll50bcZs+ccUTEfcCmERaZClyXln0IGCfpP7Wjbb2mgWMJ\nviy6YRHxy4h4NE3/FljJmzNH/n02oMFjCf59NiwiXk2Te5GNa1ePT1SHsE8bbZs9UzgaUB02dGBw\n95yUTm9/mm7/Yg2Q9C6ys7mHqj7y77OgEY4l+PfZMEljJD0C/BK4KyKWVi2ySwgbeCmFsOvqxgCg\ndd5y4MiIeFXSJ4BbgPd1uE1dT9K+ZP9i+7P0r2Vr0ijH0r/PAiJiO3CcpP2BWyQdHRFPjbDKqGdz\n/XTG0XBg0EYWEb+tnN5GxB3AW0b7F8igk7QH2V9010fErTUW8e+zQaMdS/8+mxMRvwEWAVOqPqqE\nsMmHsEfaVq8VDlG/Gi4AzgOQdBLwUkT8ql0N60F1j2W+713SRLLLtkf8IRn/C3gqIr5d53P/Phs3\n4rH077Nxkg6SNC5N7w2cDqyqWqwSwoZdQ9h19UxXlaQfAUPA2yU9A8wG9gQiIq6OiH+V9ElJTwOv\nAJ/rXGu722jHEjhb0kVkwcvXyK60sDokTQLOBZ5IfckB/DVwJP59FtLIscS/zyIOAX6QHnExBrgh\n/RZHDWGPxAFAMzMrpNe6qszMrMNcOMzMrBAXDjMzK8SFw8zMCnHhMDOzQlw4zMysEBcOa4qkcela\n+sr7yZJu62SbamlXu1LQ6kFJy1MWIf/Z1ZI+MMr6U0dbplUkrWtF0lrSXEmnjrLMZEm/v7vfZd3F\nhcOadQDw+ap53RoKarpdKTjViI+R3ab++IhYvMuXR8yKiOq0brWzgGOaaWO1dNuIkbTkv1NEzI6I\n0VLGQ8AftOL7rHu4cFizvgYclR6kMy/N20/STyStlHR9ZUFJEyQNS1oq6Y5atxOXdK2kb0taLOlp\nSdPS/F3OGCR9R1Ll1h3rJH013SV1iaTjJN2ZHkgzK7f5cZJul7RK0lW5bZ0u6X5JyyTdIGmf3Ha/\nLmkZcHZVO4+UdLekxyTdpezBQx8C5gFT0/HYq2qdRZImpOnNkr6s7GFO90s6OP2L/EzgG2n9d0s6\nKh2rpZLukfS+tP5Ryh6285ikv5O0OXec7pV0K7Aizbs5rf+EpP+ab1Kt/6Cpbd+S9GTat7en+cem\n73xU0o25W1hcm/vvtE7SnHTG9Zik90k6ErgQ+PO0X5MknZ3a84ik4VrtsB7Q6YeM+NWbL7JbQDye\nez+Z7Bkfh5D9xXQ/2b809wAWA29Py32G7Nnz1du7lux2CAAfBNbmtrsgt9x3gPPS9DpgVpr+FvAo\nsA9wEPDL3PqvpvYKWAhMA94O3APsnZb7b8Df5Lb7xTr7vQCYkaY/B9ycps8HrqizziJgQpreDnwy\nTc8D/jq3/9Ny6/wf4D1peiJwd5q+DfhMmv4T4De5/dwMHJHbxtvSn78DPAEckNu/A2u0czswPU3/\nj8r+AI8BJ6fpucC3qtuctvn5NH0RcHWang38Re47HgcOSdP7d/p37Fdzr565V5X1hCURsRFA0qPA\nu4CXgd8F7pIksrPc5+usfwtARKyU9I4Gv7NyNvIE8NbI7pr6qqTXld1GutKu9ald/0z2BMQ3gKOB\nxaldbyErdhU31Pm+3wc+naavJ/vLv4g3IuJf0/Rysi6uXUh6K1nR/UlqG6l9le+fmqZ/BPx9btUl\nEfFM7v2fSzorTR8GvBdYMkLbtgE/TtM/BG5Mx3BcZA//guyBPz+utTJwc26/Pl1nmfvI7p30Y+Cm\nEdpiXcyFw1rpjdz0NrLfl4AnI2JS7VXqrl/5C/M/2LVL9XfqrLO9av3t7Px9V/fpR9r+wog4t05b\nXqkzf3fHB7bmpivHqNoYYFNETBjl+6u7nHa0WdJk4FTgxIh4Q9Ii3nzsRlP5rkaftlc5/vX2i4j4\nvKQTgE8ByyVNiIjRnkZpXcZjHNaszcB+DSy3GjhY2a3EkbSHGntiW+Uvq/XA0ZLeIultNPBYy6r1\nAU5MYxNjyO6keh/Zs5cnSXpPatc+kt7bwHbvB85J0zOA/9tge2q1K28zsD9ARGwG1knaMb4i6ffS\n5IPsHHcZ6S6m48iKzxvKrtY6qYG2jc1t+1zgvsie4fCidl4p9sdkXXyN2rFfkI3RRMTSiJgNvMCu\nzyixHuHCYU2J7PkHiyU9rp2D47sskpbbSvaX0bzUffUIWXdLzeVrrP8sWdfIk8B84OER1qm3vSXA\nlWSDxj+PiJsj4t+BmcA/S3qMrCC8v4Ht/inwubQv5wJ/NsKytdpSb9vzgS+lweV3p21fkAaknyQb\nPAf4AvAX6fvfQ9YVWMudZA84WgF8FXiggTa8AkyU9ATZ1VB/m+afD3wzfeeHcvMb2a/bgE9XBseB\nv0+/mceBxRHxeJ31rIv5tupmPUTS3hHxWpr+LNlgdr3xhKLb3hwRjZxF2oDzGIdZbzle0pVkXV6b\ngP/Swm37X5HWEJ9xmJlZIR7jMDOzQlw4zMysEBcOMzMrxIXDzMwKceEwM7NCXDjMzKyQ/w8u32Qy\nn7ZEHwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efbe4c74ba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "class Example3():\n",
    "    def __init__(self, h=0.05):\n",
    "        self.h = h\n",
    "        self.q = 2\n",
    "        self.points = np.array([\n",
    "                [0.0, 0.0],\n",
    "                [2.0, 0.0],\n",
    "                [4.0, 1.0],\n",
    "                [4.0, 4.0],\n",
    "                [3.0, 5.0],\n",
    "                [1.0, 4.0]\n",
    "            ])\n",
    "    def __call__(self, p):\n",
    "        x = p[:, 0]\n",
    "        y = p[:, 1]\n",
    "        return x*y\n",
    "\n",
    "f = Example3()\n",
    "fe = exact_quad(f)\n",
    "print(fe)\n",
    "maxit = 3\n",
    "error = np.zeros(maxit, dtype=np.float)\n",
    "for n in range(1, maxit+1):\n",
    "    fn = polygon_quad(f, n)\n",
    "    error[n-1] = np.abs(fn-fe)\n",
    "print(error)\n",
    "print(error[:-1]/error[1:])\n",
    "\n",
    "g = plt.figure()\n",
    "axes = g.gca()\n",
    "axes.set_xlim([-0.1,4.1])  \n",
    "axes.set_ylim([-0.1,5.1])\n",
    "axes.set_axis_off()\n",
    "coord = [[0,0], [2,0], [4,1], [4,4],[3,5],[1,4]]\n",
    "coord.append(coord[0])\n",
    "xs, ys = zip(*coord)\n",
    "plt.plot(xs, ys, linewidth=1.0)\n",
    "plt.title('the figure of Integral area')\n",
    "plt.show()\n",
    "\n",
    "x = np.arange(1,maxit+1)\n",
    "y = error\n",
    "plt.plot(x,y)\n",
    "xlabel('the number of integral points')\n",
    "ylabel('error')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "35.0\n",
      "33.75\n",
      "35.0\n",
      "35.0\n",
      "[  1.25000000e+00   2.48689958e-13   2.55795385e-13]\n",
      "[  5.02633887e+12   9.72222222e-01]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEKCAYAAAACS67iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFhZJREFUeJzt3XuwrXVdx/HP90Bc5TIMDCDCIVKoUDNHMftDg0GrIbzh\nbYrk0qgVOE0UqFAjYhpdvJXVADGkfwiagzAkJqU0KCBlkY5KYxoKIndUEOUi59cfz2/B9/zOWnvd\nnsvv93ver5k97L0uz/PstQ7f9TvPe699LIQgAED5Ng19AACAdjDQAaASDHQAqAQDHQAqwUAHgEow\n0AGgEgz0kTCzzWa2xcxaec7N7HfM7A4zu9/M9jKzB8zs4Da23SUzO9TMbjSz75vZqUMfzzrM7AQz\n++zQx4F8MNArZWY3m9lRycWtvOnAzLaX9G5JR4cQdg8h3BdC2C2E8M02tt+xMyR9JoSwRwjhA+mV\nZna1mZ28yIbM7CIzO6f1I1wObyTB4xjoWMV+knaUdFPXOzKz7Vre5GZJX2l5m63r4Puetz9mQQV4\nEitkZh+SdJCkK+IpkT+cXCXpeDP7lpndZWZnuvuYmb3FzL5uZneb2SVmtueUbT9N0v/EL79rZv8a\nL99iZofEz/cysyviaY0bzOwdk1MD0079+FVxPI3wOTN7j5ndI+lt8fKTzeyrZnavmX3SzA7a4Pt/\niZl92czuM7PPmNlh8fJPSzpS0t/Ex+Wpcx7HF5rZrWZ2mpndaWa3mdmJ8brXS/oNSWfEbV0eL9/f\nzD4WH99vmNmb3PZ2MrMPxuP6ipmdbma3uutvNrMzzOyLkn5gZpvM7M3xObk/fk8v2+iYk+P/qJnd\nbmbfNbN/M7OfddddZGZ/a2afMLMHJP2Sme1gZn8Z/3zcHq/fMd5+z/ic3hWfgyvM7MmLHgt6EkLg\no8IPSTdLOtJ9vVnSFknnSdpB0jMlPSTpsHj970m6TtL+kn5C0t9J+vCMbW+W9Jgkc5c9JumQ+Pkl\nkj6sZhX/M5JukXRNct9N7r5XSzo5fn6CpEcl/a6aBceOkl4q6WuSDo2XnSnp2hnHdqikH0g6StJ2\nkk6X9L+Stk/3NeP+/lheGI/lbXFbvyrpQUl7xOsvknSOu69J+oKks+LtD5b0dUkvitefG7e/u6Qn\nS/qipFuS5+y/4nU7xsuOk7Rv/PxV8Xvb1z1W12zwvZwoaZf4fL5H0o3uuoskfVfSL8Svd5T0XkmX\nSdpD0q6SLpf0znj9XpJeHm+3q6SPSLp06D/nfCTP+dAHwEdHT2wzHI5yX08G6f7ushskvTp+/lVt\n/QKwv6RH/OCdsi0/lLdIOiQO3EckPdVd9w4tN9C/mezvSkknua83xcF64JRj+yNJl7ivTdK3Jb0g\n3deMxy0d6A8mx3qnpCPi5+lAP2LKsb9F0oXx82+o6Q6T635rykA/Yc7zeqOkY91jNXOgJ/fbMz5H\nu7lj/4fkNj+Q9JPu6+dL+r8Z23uWpHuH/nPOx9Yf2wtjc6f7/IeSnhQ/3yzp42a2JX5talan+0q6\nfYnt76Nmdfptd9mtM247S3r7zZLeb2bvdscWJB0w5bZPlvStyRchhBBPaxyw5DFM3BtC2OK+9o9Z\narOkA8zsPnecmyRd445t3uPir5eZvU7S76tZ7UvN6njveQcdT2m9S9Ir4+1D/Nhb0gPp/s1sHzWr\n+f80s8nFm+L3IDPbWdL7JP2ymhcHk/QkM7MQJzyGx0Cv17L/k92iZmV6/Zr7vVvSjyU9Rc3pBkk6\n0F3/YPzvLmpWhFITWb302G+R9CchhIsX2P93JD09uexAJYOyJelx3qpmRXvYjNt/R83jMmkQ0zrA\n49uMneB8NX9zuj5edqPikJ3j1yUdq+ZvabeY2R5qTrH4+/rjv0fNi9XhIYRpL+B/IOlpkp4bQrjb\nzH5OzemhyYsrMkAUrdcdak6BeBsNgvMkvWsSG81sHzN7yQa3n7qtuJq9VNLZZrazmf20pNe56++R\ndJuaOLspxtCfmvO9nCfpzEnUM7M9zOyVM277UUnHmNmRZrZ9DMIPSVr3hWqaO7X1Y/zvkh6IYXMn\nM9vOzA43s+fE6/9R0ltjYDxA0ilztr+rmtMk98TH6iRt+2I1y26SHlYTrneV9KfaYPDGVfYFkt4X\nV+syswPM7MVuez+SdL+Z7SXp7AWPAz1ioNfrXEl/HH+i4rR4Wfo/tP/6/Woi2FVm9n01gfSIDba/\n0bbepOav5bdL+qCaQPqwu/71an4e/B410fTajb6REMJl8fu5xMy+J+lLkn5lxm2/Jul4SR9Q87eF\nY9Scc/7xjOOe931tdP2Fkg6Pj/Gl8cXs19ScX75Z0l1qhuTu8fbnqHkxu1nSVWoGvH9cttp3COEm\nNT/v/3k1L9CHS/rcnOOb+JCav9ncJunLap7Ped6s5m9Vn4+P81VqIrPUnG7ZRc1zdp2aroHMGKe/\n0DUzO1fNT2acNPSx5MTMflvSa0IIRw59LKgDK3S0zswOM7NnxM+PUPPTHJcOe1TDM7P9zOwXrXGY\nmvPSo39c0B6iKLqwm6SLzWx/NeeZ/yKEcMXAx5SDHdT0gIMlfU/SxWp+3h9oBadcAKASnHIBgEow\n0AEgQ2baz0zXm+n+Re/DOXQAyIyZnqXmx4gPknTanJs/cT/OoQNAPsz0MjXvX7hazfs0nh2CHl3k\nvpxyAYAMmMnM9FY1b4p7laTnSTpl0WEuMdABYHBm2knNu3uPUzPIXyzpsyE8/ovdFsI5dAAYkJn2\nk/RxNb/c7QVqfpncGyQ9Y9ltsUIHgIHE+HmDpE9Jeq2aX4D2V5LeGcJSv7ZaEit0ABiEi5+nhqCP\nxMuOU/N787f5B8wX2iY/5QIA/TGTqfmXrE6R9PIQ9B/x8l3V/Mthv7nsufMJVugA0JMYPy9Q8+OI\nzwtBt7mrz9IKIdRjoANAD9L4GYJ+6K47TCuGUI8oCgAdS+NnMsxNa4RQjxU6AHRoWvxMvEJrhNCt\n9kUUBYD2zYqfyW3WDqEeK3QAaNmc+OmtHUI9BjoAtGij+JncrpUQ6hFFAaAlG8XP5HathVCPFToA\ntGCB+Om1FkK3OgaiKACsbpH4mdy+1RDqsUIHgBUtET+9VkOox0AHgBUsGj+T+7QeQj2iKAAsadH4\nmdynkxDqsUIHgCUsGT+9TkKoRxQFgAUsGz+T+3YWQj1W6AAwx4rx0+sshHoMdADYwCrxM7l/pyHU\nI4oCwAyrxM/k/p2HUI8VOgBMsUb89DoPoR5RFACcdeJnsp1eQqjHCh0Aohbip9dLCPUY6ACg9eNn\nsq3eQqhHFAUweuvGz2RbvYZQjxU6gFFrKX56vYZQjygKYJTaip/JNnsPoR4rdACj03L89HoPoR4D\nHcCotBk/k+0OEkI9oiiA0WgzfibbHSyEeqzQAYxCB/HTGyyEekRRAFXrIn4m2x80hHqs0AFUq8P4\n6Q0aQj0GOoAqdRU/k30MHkI9oiiA6nQVP5N9ZBFCPVboAKrScfz0sgihHlEUQBW6jp/JvrIJoR4r\ndADF6yl+etmEUI+BDqBofcTPZH9ZhVCPKAqgWH3Ez2R/2YVQjxU6gCL1GD+97EKoRxQFUJQ+42ey\n3yxDqMcKHUAxBoifXpYh1GOgAyhC3/Ez2Xe2IdQjigLIXt/xM9l31iHUY4UOIGsDxU8v6xDqEUUB\nZGmo+JkcQ/Yh1Ot1hW6mHULQI33uE0B5Bo6fXvYh1Ov7HPo3zfTcnvcJoCAxfl4taUc18XOQYe5C\n6OlD7H8VfQ/0xyRdY6bX9LxfAAUYMn4mx1FMCPX6HugXSvqCpD8309lm/JQNgEaMn/8i6YwQdHYI\n2jLg4RQTQr1eo6iZDpT035KeI+nDan6e9MShXoUBDC+H+JkcT1Eh1Ot1hRyCbpV0naQXSjpS0sNq\nTsEc0OdxAMhDjJ8fknScmvg56DCPigqh3hCnPM6X9MYQ9JCk10n6mKQbiKXAuOQSP70SQ6g3xED/\npKSnmOmZISiEoHMlnSrpSmIpMA65xE+v1BDq9T7QQ9CP1cTRN7jLLpN0tIilQPUyi59ekSHUG+Sd\noi6OHuhfmZNfvkMsBSqSW/z0Sg6h3iArYRdHX51cfoeIpUB1Mo2fXrEh1Bvy1Mb5kt6YXkgsBeqS\nY/z0Sg+h3pAD/fE4ml5BLAXqkGP89GoIod5gA31aHJ1yG2IpUKiM46dXfAj1Bv31ubPi6JTbEUuB\nQuQcP71aQqg36Ip3VhydcjtiKVCAAuKnV0UI9XI4hTE1jqaIpUDeco+fXk0h1MthoM+MoyliKZCn\n3OOnV1sI9QYf6IvE0Sn3IZYCmSgkfnpVhVAvi39TdNE4OuV+xFJgIKXET6/GEOplsbJdNI5OuR+x\nFBhAYfHTqy6EelkM9GihOJoilgL9Kil+erWGUC+ngb5wHE0RS4F+lBQ/vZpDqJfNQF8ljk7ZBrEU\n6EiB8dOrNoR6WUTRiVXj6JTtEEuBlpQYP70YQm+SdHyt584nslrBrhpHp2yHWAq0oOD46Z0l6Zra\nh7mU2UCPVoqjKWIpsJ5S46c3hhDq5TjQV46jKWIpsJpS46c3lhDqZTfQ24ijU7bpY+nbiaXAbIXH\nT28UIdTLKopOtBVHp2x3XzWx9DZJJ5S46gC6Unr89MYUQr0sV6ptxdEp271T0lGSfiRiKfC4SuKn\nN5oQ6mU50KNW4mgqxtITRCwFJNURP72xhVAv54HeWhxNEUuBRg3x0xtjCPWyHehdxNEp+yCWYrQq\nip/e6EKol2UUnegqjk7ZD7EUo1FT/PTGGkK9rFekXcXRKfshlmIUKoyf3ihDqJf1QI86iaMpYilq\nV1v89MYcQr0SBnpncTRFLEWtaouf3thDqJf9QO8jjk7ZJ7EU1ag0fnqjDqFe1lF0oq84OmW/xFIU\nq9b46RFCt1bEyrOvODplv8RSFKny+OmNPoR6RQz0qJc4miKWojQ1x0+PELqtkgZ6b3E0RSxFKWqO\nnx4hdLpiBvoQcXTKMRBLka0RxE+PEDpFEVF0Yqg4OuU4iKXIxhjip0cIna2oFeZQcXTKcRBLkYUR\nxU+PEDpDUQM9GiSOpoilGNpY4qdHCN1YiQN9sDiaIpZiKGOJnx4hdL7iBnoOcTTlYumfEUvRNRc/\nTx9B/PQIoXMUFUUncomjKWIpujS2+OkRQhdT5EoylziaIpaiKyONnx4hdAFFDvQoiziaIpaibWOM\nnx4hdHElD/Rs4miKWIq2jDF+eoTQ5RQ70HOMoyliKdYx4vjpEUKXUGQUncg1jqaIpVjGmOOnRwhd\nXtErxlzjaIpYikURP7dCCF1S0QM9yjKOpoilmGfs8dMjhK6mhoGebRxNEUsxy9jjp0cIXV3xA72E\nOJri1/DCI35ugxC6oqKj6EQpcTRFLB034ue2CKHrqWJlWEocTRFLx4v4ORMhdA1VDPSoiDiaIpaO\nD/FzOkLo+moa6MXE0RSxdDyIn9MRQttRzUAvMY6miKV1I35uiBDagiqi6ESpcTRFLK0L8XNjhND2\nVLUCLDWOpoil9SB+LoQQ2pKqBnpUZBxNEUvLR/ycjxDarhoHerFxNEUsLRfxcz5CaPuqG+g1xNEU\nsbQsxM+FEUJbVlUUnagljqaIpXkjfi6OENqNKld6tcTRFLE0X8TPpRFCO1DlQI+qiKMpYml+iJ/L\nIYR2p+aBXk0cTblYeoqkTxBLh0P8XA4htFvVDvQa42gqBF0u6UXi3ywdBPFzJYTQDlUZRSdqjaMp\nYmm/iJ+rIYR2r+oVXa1xNEUs7Q/xcy2E0I5VPdCjKuNoiljaPeLn6gih/RjDQK82jqaIpd0hfq6O\nENqf6gf6GOJoiljaLuLn2gihPak6ik6MJY6miKXrIX6ujxDar1Gs3MYSR1PE0tURP1tDCO3RKAZ6\nNIo4miKWLo/42Q5CaP/GNNBHE0dTxNLFufj5z5Jew2mq1RBChzGagT7GOJoilm4siZ9vD0H1B6bu\nEEIHMIooOjHWOJoilm6N+NkuQuhwRrVCG2scTRFLn0D87AQhdCCjGujRKONoilhK/OwCIXRYYxzo\no42jqSSWXmmm1w59TH0hfraPEDq80Q104ui2Yiw9WtK5Y4ilxM/OEEIHNqooOkEcna72WEr87A4h\nNA9Vr8RmIY5OV3MsJX52jhCagVEO9Ig4OkWNsZT42S1CaD7GPNCJozPUFEuJn90ihOZltAOdODpf\n6bGU+NkLQmhGRhlFJ4ijiyktlhI/+0EIzU9RK662EUcXU1IsJX72ihCamVEP9Ig4uoASYinxsz+E\n0Dwx0ImjC8s5lhI/+xNPaf21CKHZGf1AJ44uL7dYSvzs3Ssk7S9CaHZGHUUniKOrGTqWEj/7RwjN\n2+hX6BJxdFVDxlLi52AIoRljoD+BOLqCIWIp8XMYhND8MdCfQBxdUZ+xlPg5DEJoGRjoEXF0fV3H\nUuLnoAihBSCKOsTRdrQdS4mfwyKEloMVukMcbUebsZT4mQVCaCEY6NsijragjVhK/BweIbQsDPRt\nEUdbsk4sJX4OjxBaHgZ6gjjavmVjKfEzG4TQwhBFpyCOdmNeLCV+5oMQWiZW6FMQR7uxUSwlfmaH\nEFogBvpsxNEOTIulxM+8EELLxSmXGcy0vaSbJR0Tgr409PHUyEwvlfT3kh6WdIGkczhfPqx42utT\nkj4Zgt479PFgOQz0DZjpbEl7h6BThz6WWpnp6ZIOCkFXDn0skMx0nKSzJT07BD068OFgSQz0DRBH\nMSaE0PJxDn0DxFGMDCG0cKzQ5zDTsZLODEHPH/pYgK7EEHqtpGfwJqJysUKfj3eOomq8I7QeDPQ5\neOcoRoB3hFaCUy4LII6iVoTQurBCXwBxFBUjhFaEFfqCiKOoDSG0PqzQF0ccRTUIoXVioC+IOIrK\nEEIrxCmXJRBHUQNCaL1YoS+BOIpKEEIrxQp9ScRRlIwQWjdW6MsjjqJIhND6MdCXRBxFwQihleOU\nywqIoygNIXQcWKGvgDiKAhFCR4AV+oqIoygFIXQ8WKGvjjiK7BFCx4WBviLiKApBCB0RTrmsgTiK\nnBFCx4cV+hqIo8gcIXRkWKGviTiKHBFCx4kV+vqIo8gKIXS8GOhrIo4iQ4TQkeKUSwuIo8gFIXTc\nWKG3gDiKjBBCR4wVekuIoxgaIRSs0NtDHMVgCKGQGOitIY5iYIRQcMqlTcRRDIEQiglW6C0ijmIg\nhFBIYoXeOuIo+kQIhccKvX3EUfSCEIoUA71lxFH0iBCKrXDKpQPEUXSNEIppWKF3gDiKHhBCsQ1W\n6B0hjqIrhFDMwgq9O8RRtI4Qio0w0DtCHEVHCKGYiVMuHSKOok2EUMzDCr1DxFG0jBCKDbFC7xhx\nFG0ghGIRrNC7RxzFWgihWBQDvWPEUbSAEIqFcMqlB8RRrIoQimWwQu8BcRRrIIRiYazQe0IcxbII\noVgWK/T+EEexMEIoVsFA7wlxFEsihGJpnHLpEXEUiyCEYlWs0HtEHMWCCKFYCSv0nhFHsRFCKNbB\nCr1/xFFMRQjFulihD8BMn5Z0lKR/GvpYkJWdJe0n6edD0KNDHwzKw0AHgEpwygUAKsFAB4BKMNAB\noBIMdACoBAMdACrBQAeASjDQAaASDHQAqAQDHQAqwUAHgEow0AGgEgx0AKgEAx0AKsFAB4BKMNAB\noBIMdACoBAMdACrBQAeASjDQAaASDHQAqMT/A/pcgeEPXmOJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efbe4c8b7f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEPCAYAAABY9lNGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+cVPV97/HXexGxGLRosFWTaAhoGmKihIp70ThIKmsx\nq6QYRWiirY/khmhbA/ca+ogVbyu3knsTQr1tgrEWaBGSSCISkaTglqqgRjAaq2hSg0pSeTSaIiYQ\nAp/7x/csjuvs7uwyZ8/M7vv5eMwjZ2a+58xnjhM++/2tiMDMzKyjpqIDMDOz+uQEYWZmFTlBmJlZ\nRU4QZmZWkROEmZlV5ARhZmYV5ZogJN0m6SVJj3dT7ncl7ZP0kTzjMTOz6uVdg7gdmNxVAUlNwF8D\n63KOxczMeiDXBBER9wOvdFPsGuAbwM48YzEzs54ptA9C0gnAxRHxd4CKjMXMzN6o6E7qhcB1Zc+d\nJMzM6sRhBX/+OGCFJAFvBS6QtC8iVncsKMmLRpmZ9UJE9OqP776oQYhOagYRMTJ7vJPUDzGrUnIo\nK+9HjR433HBD4TH0p4fvp+9lvT4ORa41CEnLgRJwrKTngRuAw4GIiMUdiruGYGZWR3JNEBFxeQ/K\n/lGesZiZWc8U3UltBSmVSkWH0K/4ftaO72X90KG2UfUVSdEosZqZ1QtJRB13UpuZWQNygjAzs4qc\nIMzMrCInCDMzq8gJwszMKnKCMDOzipwgzMysIicIMzOryAnCzMwqcoIwM7OKnCDMzKwiJwgzM6vI\nCcLMzCpygjAzs4qcIMzMrCInCDMzq8gJwszMKmqoBPHaa0VHYGY2cDRUgvjYx+DAgaKjMDMbGBoq\nQezcCZ/7XNFRmJkNDA2VIFatghUrYOnSoiMxM+v/ck0Qkm6T9JKkxzt5/3JJ388e90s6ravrjRgB\nd98Nc+bA/ffnE7OZmSV51yBuByZ38f6/Ax+MiPcDfwXc2t0Fx4xJNYhLLoHnnqtRlGZm9ia5JoiI\nuB94pYv3N0fEf2VPNwMnVnPdlhaYOxc+/GHYtasGgZqZ2ZvUUx/EVcDaagtfcw2ccw5cdhn8+tc5\nRmVmNkAdVnQAAJImAlcCZ3dVbt68eQePS6USixaVuOCC1CexcGG+MZqZNYK2tjba2tpqci1FRE0u\n1OkHSCcBd0fE+zp5/33AnUBLRPyoi+tEpVhfeQXOOgs+8xn45CdrFbWZWf8giYhQb87tixqEsseb\n35DeQUoOf9hVcujK8OGwZg2cfTaMGgWTJh1CpGZmdlCuNQhJy4EScCzwEnADcDgQEbFY0q3AR4Dt\npCSyLyLO7ORaFWsQ7dra4NJLYeNGOPXUmn4NM7OGdSg1iNybmGqluwQB8NWvwoIFsHkzHHNMHwVm\nZlbHnCDKzJ4NW7fCunUweHAfBGZmVsecIMrs3w8XXwzHHw9f+QqoV7fFzKx/OJQEUU/zIGpi0CBY\nvjw1M3noq5lZ79XFPIhaGzYsrdnU3AynnAJTphQdkZlZ4+l3TUzlNm2C1lbYsAFO63IZQDOz/slN\nTJ1obk7NTK2taS8JMzOrXr9OEAAzZsDMmanjes+eoqMxM2sc/bqJqd2BA2lRv8MPh2XLPLLJzAYO\nNzF1o6kJ/uEfYNs2mD+/6GjMzBpDvxzFVMnQobB6NYwfn5bimDat6IjMzOrbgGhiKrd1K5x/Pqxd\nC+PG1SAwM7M65iamHjjjDFi8OHVav/hi0dGYmdWvAdPEVG7q1NQf0doK//qvcOSRRUdkZlZ/BlwT\nU7sIuOIK2L0bvv711JFtZtbfuImpF6TU1LRzJ3zuc0VHY2ZWfwZsggAYMgRWrYIVK2Dp0qKjMTOr\nLwOyD6LciBFpYb+JE2HkyLR1qZmZDfAaRLsxY1IN4pJL4Lnnio7GzKw+OEFkWlpg7lz48Idh166i\nozEzK96AHcVUSQTMmgXbt6dZ14cN+AY4M2t0HsVUIxIsWgS/+hXMmVN0NGZmxXKC6GDw4DQvYu3a\ntKe1mdlA5UaUCoYPhzVr0oimUaNg0qSiIzIz63u51iAk3SbpJUmPd1FmkaRnJT0m6fQ84+mJ0aNh\n5Uq4/PK0LIeZ2UCTdxPT7cDkzt6UdAHwrogYDXwS+HLO8fRIqQQ33ZRGNr38ctHRmJn1rVwTRETc\nD7zSRZGLgKVZ2YeAoyX9Vp4x9dRVV6UEMW0a7NtXdDRmZn2n6E7qE4EXyp7vyF6rKwsWpBVfP/3p\nNBTWzGwgaKhO6nnz5h08LpVKlEqlPvncQYNg+XKYMAEWLoRrr+2TjzUz67G2tjba2tpqcq3cJ8pJ\nOgm4OyLeV+G9LwP3RcTK7PnTwLkR8VKFsrlPlOvO9u3Q3Ay33gpTphQaiplZVep9opyyRyWrgY8B\nSDoL+Hml5FAvTjoJ7rwz7SPxxBNFR2Nmlq9caxCSlgMl4FjgJeAG4HAgImJxVuYWoAV4DbgyIrZ0\ncq3CaxDt/umf0h4SDz0Exx1XdDRmZp07lBqE12Lqpeuvh/XrYcMGOOKIoqMxM6vMCaIABw7ApZem\nTYeWLUvrOJmZ1Zt674Pol5qaYMmSNMt6/vyiozEzq72GGuZab4YOTcuCjx8Pp56aJtOZmfUXbmKq\nga1b4fzz0wqw48YVHY2Z2evcxFSwM86AxYvh4othx46iozEzqw03MdXI1KmpP6K1FTZuTEtzmJk1\nMjcx1VBEmkS3e3fadKjJ9TMzK5ibmOqElJqadu5ME+nMzBqZE0SNDRkCq1bBihWwdGnR0ZiZ9Z77\nIHIwYgTcfTdMnAgjR6atS83MGo1rEDkZMybVIC65BJ57ruhozMx6zgkiRy0tMHdu2pFu166iozEz\n6xmPYspZBMyalfaSWL0aDnOjnpn1IY9iqmMSLFoEv/oVzJlTdDRmZtVzgugDgweneRFr18JXvlJ0\nNGZm1XGDRx8ZPhzWrEkjmkaNgkmTio7IzKxrrkH0odGjYeVKuPzytCyHmVk9c4LoY6US3HRTGtn0\n8stFR2Nm1jmPYirI7NlpmfB161IfhZlZHrzlaAPavz8tD3788anj2luWmlkePMy1AQ0aBMuXw+bN\nsHBh0dGYmb2ZRzEVaNiwtGZTczOccgpMmVJ0RGZmr3MTUx3YtCltNLRhA5x2WtHRmFl/UtdNTJJa\nJD0t6RlJ11V4/+2SNkjaIukxSRfkHVO9aW5OzUytrWkvCTOzepBrDUJSE/AMMAn4CfAIcFlEPF1W\n5ivAloj4iqTfAe6JiHdWuFa/rUG0u/56WL8+1SSOOKLoaMysP6jnGsSZwLMRsT0i9gErgIs6lDkA\nHJUd/yawI+eY6taNN8KJJ8JVV6VF/szMipR3gjgReKHs+YvZa+VuBP5Q0gvAGuCanGOqW01NsGRJ\nmmU9f37R0ZjZQFcPo5imA7dHxBclnQX8IzCmUsF58+YdPC6VSpRKpb6Ir08NHZqWBR8/Hk49FaZN\nKzoiM2skbW1ttLW11eRaefdBnAXMi4iW7PlngYiIm8vK/ACYHBE7suc/AsZHxH92uFa/74Mot3Ur\nnH9+WgF23LiiozGzRlXPfRCPAKMknSTpcOAyYHWHMtuBDwFkndRDOiaHgeiMM2Dx4jTbeseA7ZUx\nsyLl2sQUEfslXQ18h5SMbouIpyTdCDwSEWuAOcCtkq4ldVh/PM+YGsnUqak/orUVNm6EI48sOiIz\nG0g8Ua7ORcAVV8Du3WnToSYvjmJmPVDPTUx2iKTU1LRzZ5onYWbWV7pNEJIGZc0/VpAhQ2DVKrjj\nDli2rOhozGygqKqJSdLDEXFmH8TTVQwDsomp3JNPwsSJKVmcfXbR0ZhZI8h9PwhJXwQGAyuB19pf\nj4gtvfnQ3nCCSO69F668Eh58EN75pgVJzMzeqC8SxH0VXo6IOK83H9obThCvW7Qo9Us8+CAcdVT3\n5c1s4PKOcgNMBMyaBdu3p1nXh9XDfHgzq0u5j2KSdLSkL0j6Xvb4v5KO7s0H2qGTUi3iV7+COXOK\njsbM+qtqh7n+PfAq8NHssQu4Pa+grHuDB6d5EWvXpj2tzcxqrdo+iMci4vTuXsuTm5gqe/bZNKJp\n+XKYNKnoaMys3vTFRLlfSjo4sFLSBOCXvflAq63Ro2HlSrj88rQsh5lZrVRbg3g/sBRo73d4Bfh4\nRDyeY2wdY3ANogtf/SosWACbN8MxxxQdjZnVi1xHMWXbhk6LiK9JOgogInb15sMOhRNE92bPTsuE\nr1uX+ijMzPpiHsT3IqLQXQmcILq3f39aHvz441PHtXr1kzCz/qQv+iD+WdIcSW+XdEz7ozcfaPkZ\nNCh1Vm/eDAsXFh2NmTW6amsQz1V4OSJiZO1D6jQG1yCqtH07NDfDrbfClClFR2NmReqLPojmiHig\nNx9QK04QPbNpU9poaMMGOO20oqMxs6Lk2sQUEQeAW3pzcStOc3NqZmptTXtJmJn1VLV9EOsl/YHk\nbs9GMmMGzJyZOq737Ck6GjNrNNX2QbwKDAX2A3sAkfog+mwtUTcx9c6BA3DppWnToWXLPLLJbKDp\ni1FMRwNXAH+VJYUxwO/15gOtbzU1wZIlaZb1/PlFR2NmjaTaBPH/gLOA6dnzV3G/RMMYOhTuuivN\njfjGN4qOxswaRbU7CYyPiLGStgJExCuSDs8xLquxE05ISeL88+Hkk2FcodMezawRVFuD2CdpEBAA\nkkYAB3KLynJxxhlpJ7qLL4YdO4qOxszqXbUJYhHwTeA4STcB9wNVtWhLapH0tKRnJF3XSZmPSnpS\n0hOS/rHKmKwXpk6Fq69Ow19fe6378mY2cFW95aikdwOTSCOY1kfEU1Wc0wQ8k533E+AR4LKIeLqs\nzChgJTAxInZJemtE/GeFa3kUU41EwBVXwO7dadOhpmr/TDCzhlO3e1JLOgu4ISIuyJ5/ljQ89uay\nMjcD2yLi77u5lhNEDe3dCx/6EHzwg3DTTUVHY2Z56Ythrr11IvBC2fMXs9fKnQKcKul+SQ9Kmpxz\nTEaaF7FqFdxxR5ofYWbWUbWjmPJ0GDAK+CDwDmCjpPdW2nNi3rx5B49LpRKlUqmPQuyfRoyAu++G\niRNh5EiYMKHoiMzsULW1tdHW1laTa/VFE9O8iGjJnldqYvo7YHNELMme/zNwXUQ82uFabmLKyb33\nwpVXpgX+Tj656GjMrJbquYnpEWCUpJOyeROXAas7lPkWMBFA0luB0cC/5xyXlWlpgblz4cILYVef\n7xVoZvUq1wQREfuBq4HvAE8CKyLiKUk3SrowK7MO+JmkJ4H1wJyIeCXPuOzNrrkGzjkHLrsMfv3r\noqMxs3qQaxNTLbmJKX/79sEFF8B73+sd6cz6i3puYrIGMnhwmhexdm1at8nMBrZ6GMVkdWT4cFiz\nBs4+G0aNgkmTio7IzIriGoS9yejRsHIlXH55WibczAYmJwirqFRKM6w//GF4+eWiozGzIriT2ro0\nezZs3Qrr1qU+CjNrLHW7FlMtOUEUY//+tDz48cenjmtvWWrWWDyKyXIzaBAsXw6bN3voq9lA41FM\n1q1hw9KaTc3NcMopMGVK0RGZWV9wE5NVbdOmtNHQhg1w2mlFR2Nm1XATk/WJ5ubUzNTaCjt3Fh2N\nmeXNCcJ6ZMYMmDkzdVzv2VN0NGaWJzcxWY8dOACXXpo2HVq2zCObzOqZm5isTzU1wZIlaZb1/PlF\nR2NmefEoJuuVoUPhrrtg/Hg49VSYNq3oiMys1tzEZIdkyxaYPDmtADtuXNHRmFlHbmKywowdC4sX\np07rHTuKjsbMaslNTHbIpk5N/RGtrbBxIxx5ZNERmVktuInJaiICrrgCdu9Omw41uW5qVhfcxGSF\nk1JT086dcP31RUdjZrXgBGE1M2QIrFoFd9yR5keYWWNzH4TV1IgRaWG/iRNh5EiYMKHoiMyst1yD\nsJobMwaWLk1zI37846KjMbPecoKwXLS0wNy5cOGFsGtX0dGYWW/kniAktUh6WtIzkq7rotwfSDog\naWzeMVnfuOYaOOccmD497UxnZo0l1wQhqQm4BZgMjAGmS3p3hXJvAf4E2JxnPNa3JFi0CPbuhTlz\nio7GzHoq7xrEmcCzEbE9IvYBK4CLKpT7S+Cvgb05x2N9bPDgNC/innvSntZm1jjyThAnAi+UPX8x\ne+0gSWcAb4uItTnHYgUZPhzWrIG/+AtYv77oaMysWoUOc5Uk4AvAx8tf7qz8vHnzDh6XSiVKpVJe\noVmNjR4NK1emfSQ2bkwrwJpZ7bW1tdHW1laTa+W61Iaks4B5EdGSPf8sEBFxc/b8KOCHwG5SYvht\n4GdAa0Rs6XAtL7XRD3z1q7BgAWzeDMccU3Q0Zv3foSy1kXeCGARsAyYBPwUeBqZHxFOdlL8P+ExE\nbK3wnhNEPzF7NmzdCuvWpT4KM8tP3a7FFBH7gauB7wBPAisi4ilJN0q6sNIpdNHEZP3DggVpxddP\nfzot8mdm9cmruVohXn01LcNx5ZVw7bVFR2PWfx1KDcJrMVkhhg1LazY1N8Mpp8CUKUVHZGYduQZh\nhdq0KW00tGEDnHZa0dGY9T912wdh1p3mZli4MCWJnTuLjsbMyjlBWOFmzICZM9O+1nv2FB2NmbVz\nE5PVhQMH0iS6IUPSZkPyWDazmnATkzW8piZYsgS2bYP584uOxszAo5isjgwdCnfdBePHp6U4pk0r\nOiKzgc1NTFZ3tmyByZNh7VoYN67oaMwam5uYrF8ZOxYWL06d1jt2FB2N2cDlJiarS1Onpv6I1ta0\n+uuRRxYdkdnA4yYmq1sRcMUVsHt32nSoyfVdsx5zE5P1S1Jqatq5E66/vuhozAYeJwira0OGwKpV\ncMcdaX6EmfUd90FY3RsxIi3sN3EijByZVoE1s/y5BmENYcwYWLo0zY348Y+LjsZsYHCCsIbR0gJz\n58KFF8KuXUVHY9b/eRSTNZQImDULnn8eVq+GQYOKjsisvnkUkw0YEixaBHv3wpw5RUdj1r85QVjD\nGTw4zYu45540DNbM8uFRTNaQhg+HNWvgnHNg1Cg477yiIzLrf1yDsIY1ejSsWAHTp8MzzxQdjVn/\n4wRhDa1UgptuSiObXn656GjM+hePYrJ+YfZs2LoV1q1LfRRmltT1KCZJLZKelvSMpOsqvH+tpCcl\nPSbpu5LenndM1v8sWJBWfP30p9NQWDM7dLkmCElNwC3AZGAMMF3SuzsU2wJ8ICJOB+4EPp9nTNY/\nDRoEy5fD5s2wcGHR0Zj1D3nXIM4Eno2I7RGxD1gBXFReICL+JSL2ZE83AyfmHJP1U8OGpTWbPv95\n+Pa3i47GrPHlnSBOBF4oe/4iXSeAPwbW5hqR9WsnnQR33pn2kXjiiaKjMWtsdTMPQtJM4APAuZ2V\nmTdv3sHjUqlEqVTKPS5rPM3NqZmptRUeegiOO67oiMz6TltbG21tbTW5Vq6jmCSdBcyLiJbs+WeB\niIibO5T7EPAl4IMR8bNOruVRTNYj118P69fDhg1wxBFFR2NWjEMZxZR3ghgEbAMmAT8FHgamR8RT\nZWXOAL4OTI6IH3VxLScI65EDB+DSS9OmQ8uWpXWczAaauh3mGhH7gauB7wBPAisi4ilJN0q6MCu2\nADgS+LqkrZK+lWdMNnA0NcGSJbBtG8yfX3Q0Zo3HE+Ws3/vJT2D8ePjiF9OGQ2YDSd02MdWSE4Qd\nii1bYPJkWLsWxo0rOhqzvlO3TUxm9WLs2LQ0+MUXw44dRUdj1hjqZpirWd6mTk39Ea2tsHFjWprD\nzDrnJiYbUCLSJLrdu9OmQ02uQ1s/5yYmsypJqalp5840T8LMOucEYQPOkCGwahXccUeaH2FmlbkP\nwgakESPSwn4TJ8LIkTBhQtERmdUf1yBswBozBpYuTXMjfvzjoqMxqz9OEDagtbTA3Llpy9Jdu4qO\nxqy+eBSTDXgRMGsWPP88rF6dNh8y6y88isnsEEiwaBHs3Qtz5hQdjVn9cIIwAwYPTvMi7rknDYM1\nM49iMjto+HBYswbOOQdGjYLzzis6IrNiuQZhVmb0aFixAqZPh2eeKToas2I5QZh1UCrBTTelkU0v\nv1x0NGbF8Sgms07Mng2PPQb33pv6KMwakfeDMMvB/v1pefATToAvf9lbllpj8jBXsxwMGgTLl8Om\nTfClLxUdjVnf8ygmsy4MG5bWbGpuTh3YU6YUHZFZ33ETk1kVNm1KGw1t2ACnnVZ0NGbVcxOTWc6a\nm2HhwpQkdu4sOhqzvuEEYValGTNg5szUcb1nT9HRmOUv9wQhqUXS05KekXRdhfcPl7RC0rOSNkl6\nR94xmfXWjTfCiSfCVVelRf7M+rNcE4SkJuAWYDIwBpgu6d0div0x8HJEjAYWAgvyjMmStra2okNo\nSE1NsGQJbNsG8+e//rrvZ+34XtaPvGsQZwLPRsT2iNgHrAAu6lDmImBJdvwNYFLOMRn+P+GhGDoU\n7rorzY34xjfSa76fteN7WT/yHuZ6IvBC2fMXSUmjYpmI2C/p55KOiQgvcmB164QTUpKYPBlOPrno\naMzyUY/zIDxf1RrC2LFpafALL4SjjoJHHy06ov5h2zbfy3qR6zwISWcB8yKiJXv+WSAi4uayMmuz\nMg9JGgT8NCKOq3AtdwmamfVCb+dB5F2DeAQYJekk4KfAZcD0DmXuBj4OPARcAmyodKHefkEzM+ud\nXBNE1qdwNfAdUof4bRHxlKQbgUciYg1wG7BM0rPAz0hJxMzMCtYwS22YmVnfqquZ1JJuk/SSpMe7\nKLMom1T3mKTT+zK+RtPd/ZR0bjZqbEv2+Fxfx9goJL1N0gZJT0p6QtKfdFLOv88qVHM//fusnqQh\nkh6StDW7nzdUKNPzSckRUTcP4GzgdODxTt6/APh2djwe2Fx0zPX8qOJ+ngusLjrORngAvw2cnh2/\nBdgGvLtDGf8+a3s//fvs2T0dmv3vIGAzcGaH9z8F/G12fCmwortr1lUNIiLuB17poshFwNKs7EPA\n0ZJ+qy9ia0RV3E/wsOKqRMR/RMRj2fFu4CnSHJ5y/n1Wqcr7Cf59Vi0ifpEdDiH1L3fsP+jxpOS6\nShBV6DjxbgeVf1RWvbOyaum3Jb2n6GAagaSTSTWzhzq85d9nL3RxP8G/z6pJapK0FfgP4LsR8UiH\nIm+YlAz8XNIxXV2zHifKWd95FDgpIn4h6QLgW8ApBcdU1yS9hfTX159mf/naIejmfvr32QMRcQA4\nQ9JRwLckvSci/q2LU7qtnTVaDWIH8Pay52/LXrNeiIjd7dXSiFgLDO7uL4qBTNJhpH/MlkXEXRWK\n+PfZA93dT/8+eycidgH3AS0d3nqR7PeZTUo+KrpZ0qgeE4ToPLOtBj4GB2dp/zwiXuqrwBpUp/ez\nvH1c0pmkYc9eA6tzfw/8W0R0tkO1f5890+X99O+zepLeKuno7Pg3gN8Dnu5QrH1SMnQxKblcXTUx\nSVoOlIBjJT0P3AAcTlqeY3FE3CPp9yX9EHgNuLK4aOtfd/cTmCbpU8A+4JekkQ1WgaQJwAzgiayd\nN4A/B07Cv88eq+Z+4t9nTxwPLMm2WGgCVma/x0OalOyJcmZmVlE9NjGZmVkdcIIwM7OKnCDMzKwi\nJwgzM6vICcLMzCpygjAzs4qcIKxLko7OxqK3Pz9X0t1FxlRJX8WVTUjaLOnRbCx/+XuLJb27m/Mv\n6q5MrUh6rhYzjyXdKOm8bsqcK6n5UD/L6osThHVnODCrw2v1Onmm13FlE4yq8SHS8ukfiIgH3vDh\nEZ+IiI6zVzu6GBjTmxg7ypZL6EpN/jtFxA0R0d2s2xLw32rxeVY/nCCsO/8bGJlt2HJz9towSV+X\n9JSkZe0FJY2V1CbpEUlrKy11Lel2SV+S9ICkH0r6SPb6G2oAkv5GUvuyFc9Jmp+t6vmwpDMk3Ztt\nfPKJsssfLWmNpKcl/W3ZtX5P0oOSvidppaShZdf9a0nfA6Z1iPMkSeslfV/Sd5U2uHk/cDNwUXY/\nhnQ45z5JY7PjVyX9ldLGQQ9KGpH9hd0KLMjOf6ekkdm9ekTSv0g6JTt/pNKmLt+X9JeSXi27Txsl\n3QU8mb32zez8JyRdVR5Spf+gWWxfkPSD7Lsdm71+evaZj0m6s2zphtvL/js9J2leVoP6vqRTlPac\n/+/An2Xfa4KkaVk8WyW1VYrDGkDRm1z4Ud8P0tIHj5c9P5e0x8TxpH+AHiT95XgY8ABwbFbuo6Q9\nyDte73bSMgAAvwM8W3bd1WXl/gb4WHb8HPCJ7PgLwGPAUOCtwH+Unf+LLF6R9kH/CHAs8C/Ab2Tl\n/ifwubLrzunke68GZmbHVwLfzI4/Dizq5Jz7gLHZ8QHg97Pjm4E/L/v+Hyk755+Bd2XHZwLrs+O7\ngY9mx58EdpV9z1eBd5Rd4zez/z0CeAIYXvb9jqkQ5wHgsuz4+vbvA3wfODs7vhH4QseYs2vOyo4/\nBSzOjm8APlP2GY8Dx2fHRxX9O/ajd4+6WovJGsbDEfFTAEmPAScD/wW8F/iuJJFqpz/p5PxvAUTE\nU5KOq/Iz22sXTwBHRlrl8xeS9igtb9we1/YsrjtIO+rtBd4DPJDFNZiU1Nqt7OTzmoGp2fEy0j/y\nPbE3Iu7Jjh8lNU29gaQjScn161lsZPG1f/5F2fFy4PNlpz4cEc+XPf8zSRdnx28DRgMPdxHbfuBr\n2fE/Andm9/DoSJtMQdpY5muVTga+Wfa9pnZS5n7S2kBfA1Z1EYvVMScI6429Zcf7Sb8jAT+IiAmV\nT+n0/PZ/GH/NG5s8j+jknAMdzj/A67/jjm3ukV3/OxExo5NYXuvk9UNtv99Xdtx+jzpqAl6JiLHd\nfH7HpqKDMUs6FzgPGB8ReyXdx5vvXXfaP6va3dva739n34uImCXpd4ELgUcljY2I7nY3tDrjPgjr\nzqvAsCrKbQNGKC1zjaTDVN0OYO3/KG0H3iNpsKTfpIrtEDucDzA+6ztoIq38eT9pb94Jkt6VxTVU\n0ugqrvsgMD07ngn8a5XxVIqr3KvAUQAR8SrwnKSD/R+S3pcdbub1fpGuVt08mpRk9iqNjjqritgG\nlV17BnB/pD0EXtbrI7P+kNQ0V62D3wtSH0pEPBIRNwA7eeM+GdYgnCCsS5HW339A0uN6vZP6DUWy\ncvtI/+hAkRbaAAABA0lEQVTcnDU7bSU1k1QsX+H8F0lNGj8AVgBbujins+s9DNxC6rz9UUR8MyL+\nE7gCuEPS90n/8J9axXX/BLgy+y4zgD/tomylWDq79grgf2SdvO/Mrv3HWcfwD0id2ADXAp/JPv9d\npCa8Su4lbaTzJDAf2FRFDK8BZ0p6gjT66H9lr38c+D/ZZ76/7PVqvtfdwNT2Tmrg89lv5nHggYh4\nvJPzrI55uW+zOiTpNyLil9nxpaRO5c7a+3t67VcjoppaoQ1w7oMwq08fkHQLqanqFeCPanht/1Vo\nVXENwszMKnIfhJmZVeQEYWZmFTlBmJlZRU4QZmZWkROEmZlV5ARhZmYV/X9qOJ/FqspBSgAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7efc0e6a5780>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "class Example4():\n",
    "    def __init__(self, h=0.05):\n",
    "        self.h = h\n",
    "        self.q = 2\n",
    "        self.points = np.array([\n",
    "                [1.0, 1.0],\n",
    "                [3.0, 1.0],\n",
    "                [5.0, 5.0],\n",
    "                [2.0, 2.0],\n",
    "                [0.0, 4.0]\n",
    "            ])\n",
    "    def __call__(self, p):\n",
    "        x = p[:, 0]\n",
    "        y = p[:, 1]\n",
    "        return x*y\n",
    "\n",
    "f = Example4()\n",
    "fe = exact_quad(f)\n",
    "print(fe)\n",
    "maxit = 3\n",
    "error = np.zeros(maxit, dtype=np.float)\n",
    "for n in range(1, maxit+1):\n",
    "    fn = polygon_quad(f, n)\n",
    "    error[n-1] = np.abs(fn-fe)\n",
    "    print(fn)\n",
    "print(error)\n",
    "print(error[:-1]/error[1:])\n",
    "\n",
    "g = plt.figure()\n",
    "axes = g.gca()\n",
    "axes.set_axis_off()\n",
    "coord = [[1,1], [3,1], [5,5], [2,2],[0,4]]\n",
    "coord.append(coord[0])\n",
    "xs, ys = zip(*coord)\n",
    "plt.plot(xs, ys, linewidth=1.0)\n",
    "plt.title('the figure of Integral area')\n",
    "plt.show()\n",
    "\n",
    "x = np.arange(1,maxit+1)\n",
    "y = error\n",
    "plt.plot(x,y)\n",
    "xlabel('the number of integral points')\n",
    "ylabel('error')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from IPython.display import display\n",
    "from sympy import *\n",
    "import numpy as np\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIsAAAAUBAMAAABPB9NaAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAMnZUZs0Qu91E7yKJ\nmaurDqYVAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACoUlEQVQ4EY1UT2jTUBz+krZJ1mRt3fDioYtT\nQZhbxxAUT6HgSYfxqILGwRBEtAxFcGCLHjx4sEcRYcXDQFCpOJx/DtahB9HBBhteBKvILgpTkcnU\ndX4vS5om8+CDJN/3vd/78vv9Xl6A6JgICfFqiHpEt1rVfR6RbgWqXgiwQJ99usn0ETAZQKKU9yrp\nTyAnIq8f9qe6aj4CjgSQSLFD1CWDEckwIwJpfCasTYWpYFcjklKJCKRaJOOTayH6hmaoXmpCDxzy\nnvczzZlHTbQGEmVIW/L2ZAMb73b0mRRl5iuvcsxgqO/xDkrTVDt2btOLdenwbPdHMvTyyjFmRek1\n8cSBMYMTQAlHgSum/puzqgW86ezJP4VSi1/SbGAve1iAhkQd2gieiXrOAsbzgYmeqhGrIMde2bjh\nwMRF4ALwgxFaP6R+o7qZgY6xFOOi3UDbOGSodagVpC0GXQY68VrO4EG7jfNArATtZw+EzTngOyOS\nZd5ULDAFJArE6OJ1baQKlTY20ibpbV4opBw4SVOskr5Ceck20IamTZss6+Rwl7g2p46tQKVNIbCR\n6hpDxhzpm2vzEPHlFhsWBcxhiXd0ZcT9OqA7OFBWAxsWhXYr7QDzkCssuoQB4FOLjWixsoRlFuzk\nIDOSLY6ZSFlqYMMW44yTLotuUBYt3uNgVFTkFcW9g1zCHebbX4Q4pu9oMw61yor8oqYpb0eyBrzF\nmMXMariZ31/Ori5kVz/s+lV38xOb/gU43T00xxAcpPhi6ytpsXF8sdFZvEeJdXLTUzUGTb2vAokM\nhfCIHob1n3X4MMxzefSrpjQYdoVhRgTS5tFM2Qq7iKn1ISqTbB3DrcTHox7QLNX8948i+tvyV/gO\n7nPSY23ds0TRN7uTYneC8Z8/0b8ZeqnQGLWTXgAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\sin{\\left (\\pi x \\right )} \\sin{\\left (\\pi y \\right )}$$"
      ],
      "text/plain": [
       "sin(π⋅x)⋅sin(π⋅y)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x, y = symbols('x, y')\n",
    "f = sin(pi*x)*sin(pi*y)\n",
    "f1 = f.series(x=(x, y), x0=(0, 0), n=4)\n",
    "display(f1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
